Analysis of the CBS pesticide experiment
Shihong Jia
Oct. 5, 2017
Load packages
library(ggplot2)
library(cowplot)
library(lme4)
library(pbkrtest)
library(tidyr)
library(dplyr)
library(optimx)
library(broom)Data description
Five censuses (15fa, 16sp, 16fa, 17sp and 17fa) were conducted. Information of neigboring trees (area of conspecific, heterospecific and overall species) 20 radius around the focal quadrats are added to the dataset.
Survival analysis
getwd()## [1] "E:/R/My code/Git/Changbaishan/Changbaishan1"rm(list=ls())
## Read in data
pestdat.ag_adult <- read.csv("data/pestdat.ag_adult.csv")
pestdat.ag_adult <- pestdat.ag_adult[,-c(1)]
pestdat.ag_adult$pesticide <- factor(pestdat.ag_adult$pesticide, levels=c('W', 'F', 'I', 'C'))
pestdat.ag_adult$census <- factor(pestdat.ag_adult$census)
pestdat.ag_adult$exclosure <- factor(pestdat.ag_adult$exclosure)
pestdat.ag_adult$site <- factor(pestdat.ag_adult$site)
summary(pestdat.ag_adult)## site quadrat sp. pesticide exclosure growth.form
## 1:1732 Min. :101.0 FRMA : 765 W:1087 0:2269 shrub:1741
## 2:1620 1st Qu.:202.0 TIAM : 673 F:1147 1:2063 tree :2591
## 3: 980 Median :304.0 PHSC : 327 I:1065
## Mean :304.9 EUMA : 302 C:1033
## 3rd Qu.:409.0 ACPS : 296
## Max. :510.0 ACBA : 260
## (Other):1709
## census deaths total survs
## 16fa:1169 Min. : 0.0000 Min. : 1.000 Min. : 0.000
## 16sp: 625 1st Qu.: 0.0000 1st Qu.: 1.000 1st Qu.: 1.000
## 17fa:1406 Median : 0.0000 Median : 2.000 Median : 1.000
## 17sp:1132 Mean : 0.5919 Mean : 2.795 Mean : 2.203
## 3rd Qu.: 1.0000 3rd Qu.: 3.000 3rd Qu.: 3.000
## Max. :16.0000 Max. :23.000 Max. :20.000
##
## quad.unique plot quad.sp A.sum
## 1_1_407: 37 1_0:859 1_0_101_ACKO: 4 Min. : 3990
## 1_1_503: 36 1_1:873 1_0_101_CEOR: 4 1st Qu.: 4991
## 2_1_405: 35 2_0:864 1_0_101_EUMA: 4 Median : 5452
## 1_0_510: 34 2_1:756 1_0_101_FRMA: 4 Mean : 5631
## 2_0_510: 33 3_0:546 1_0_102_CEOR: 4 3rd Qu.: 6084
## 1_1_408: 32 3_1:434 1_0_102_DEGL: 4 Max. :10411
## (Other):4125 (Other) :4308
## A.con A.het con.dens A.con_curt
## Min. : 0.000 Min. : 2374 Min. : 0.00 Min. : 0.000
## 1st Qu.: 0.000 1st Qu.: 4487 1st Qu.: 0.00 1st Qu.: 0.000
## Median : 4.235 Median : 5182 Median : 1.00 Median : 1.618
## Mean : 427.093 Mean : 5204 Mean : 5.41 Mean : 4.295
## 3rd Qu.: 658.840 3rd Qu.: 5764 3rd Qu.: 9.00 3rd Qu.: 8.701
## Max. :5308.942 Max. :10411 Max. :46.00 Max. :17.445
##
## A.het_curt A.sum_curt con.dens_curt dens
## Min. :13.34 Min. :15.86 Min. :0.000 Min. :-0.30365
## 1st Qu.:16.49 1st Qu.:17.09 1st Qu.:0.000 1st Qu.:-0.30365
## Median :17.30 Median :17.60 Median :1.000 Median :-0.13447
## Mean :17.23 Mean :17.74 Mean :1.044 Mean : 0.00000
## 3rd Qu.:17.93 3rd Qu.:18.26 3rd Qu.:2.080 3rd Qu.: 0.03472
## Max. :21.84 Max. :21.84 Max. :3.583 Max. : 3.41845
## ### Intial plots
ggplot(data=pestdat.ag_adult,
aes(x=pesticide, y=survs/total, colour=as.factor(exclosure))) +
geom_boxplot() +
facet_wrap(~growth.form)
group_by(pestdat.ag_adult, pesticide, exclosure, growth.form) %>%
summarise(surv = mean(survs/total),
se.survs=sd(survs/total)/sqrt(length(survs/total))) %>%
ggplot(aes(x=pesticide, y=surv, ymin=surv-se.survs, ymax=surv+se.survs,
colour=as.factor(exclosure))) +
labs(x="Pesticide", y="Survival") +
geom_errorbar(position=position_dodge(width = .5), width = 0.2) +
geom_point(position=position_dodge(width = .5)) +
facet_wrap(~growth.form)
###------------------------- Initial Survival models -------------------------###
pestdat.ag_adult$site <- factor(pestdat.ag_adult$site)
contrasts(pestdat.ag_adult$site) <- "contr.sum"
## Conspecific and heterospecific area of adult trees
mod.tree_adult1 <- glmer(cbind(survs, deaths) ~ site + pesticide + exclosure + census +
A.con_curt + A.het_curt +
(1|quad.unique) + (1|sp.),
data=subset(pestdat.ag_adult, growth.form=='tree'), family=binomial,
control=glmerControl(optimizer="optimx",
optCtrl=list(method=c('bobyqa', 'Nelder-Mead'))))
# sum area of adult trees
mod.tree_adult2 <- glmer(cbind(survs, deaths) ~ site + pesticide + exclosure + census + A.sum_curt +
(1|quad.unique) + (1|sp.),
data=subset(pestdat.ag_adult, growth.form=='tree'), family=binomial,
control=glmerControl(optimizer="optimx",
optCtrl=list(method=c('bobyqa', 'Nelder-Mead'))))
## conspecific density of adult trees
mod.tree_adult3 <- glmer(cbind(survs, deaths) ~ site + pesticide + exclosure + census + con.dens_curt +
(1|quad.unique) + (1|sp.),
data=subset(pestdat.ag_adult, growth.form=='tree'), family=binomial,
control=glmerControl(optimizer="optimx",
optCtrl=list(method=c('bobyqa', 'Nelder-Mead'))))
anova(mod.tree_adult1, mod.tree_adult2, mod.tree_adult3)## Data: subset(pestdat.ag_adult, growth.form == "tree")
## Models:
## mod.tree_adult2: cbind(survs, deaths) ~ site + pesticide + exclosure + census +
## mod.tree_adult2: A.sum_curt + (1 | quad.unique) + (1 | sp.)
## mod.tree_adult3: cbind(survs, deaths) ~ site + pesticide + exclosure + census +
## mod.tree_adult3: con.dens_curt + (1 | quad.unique) + (1 | sp.)
## mod.tree_adult1: cbind(survs, deaths) ~ site + pesticide + exclosure + census +
## mod.tree_adult1: A.con_curt + A.het_curt + (1 | quad.unique) + (1 | sp.)
## Df AIC BIC logLik deviance Chisq Chi Df
## mod.tree_adult2 13 5302.2 5378.4 -2638.1 5276.2
## mod.tree_adult3 13 5299.0 5375.2 -2636.5 5273.0 3.1694 0
## mod.tree_adult1 14 5269.1 5351.2 -2620.6 5241.1 31.9072 1
## Pr(>Chisq)
## mod.tree_adult2
## mod.tree_adult3 < 2.2e-16 ***
## mod.tree_adult1 1.617e-08 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1## interaction terms
mod.tree_adult4 <- glmer(cbind(survs, deaths) ~ site + pesticide*exclosure*census +
A.con_curt + A.het_curt +
(1|quad.unique) + (1|sp.),
data=subset(pestdat.ag_adult, growth.form=='tree'), family=binomial,
control=glmerControl(optimizer="optimx",
optCtrl=list(method=c('bobyqa', 'Nelder-Mead'))))
mod.tree_adult5 <- glmer(cbind(survs, deaths) ~ site + pesticide*exclosure + census +
A.con_curt + A.het_curt +
(1|quad.unique) + (1|sp.),
data=subset(pestdat.ag_adult, growth.form=='tree'), family=binomial,
control=glmerControl(optimizer="optimx",
optCtrl=list(method=c('bobyqa', 'Nelder-Mead'))))
mod.tree_adult6 <- glmer(cbind(survs, deaths) ~ site + (pesticide + exclosure)* census +
A.con_curt + A.het_curt +
(1|quad.unique) + (1|sp.),
data=subset(pestdat.ag_adult, growth.form=='tree'), family=binomial,
control=glmerControl(optimizer="optimx",
optCtrl=list(method=c('bobyqa', 'Nelder-Mead'))))
anova(mod.tree_adult4, mod.tree_adult5, mod.tree_adult6)## Data: subset(pestdat.ag_adult, growth.form == "tree")
## Models:
## mod.tree_adult5: cbind(survs, deaths) ~ site + pesticide * exclosure + census +
## mod.tree_adult5: A.con_curt + A.het_curt + (1 | quad.unique) + (1 | sp.)
## mod.tree_adult6: cbind(survs, deaths) ~ site + (pesticide + exclosure) * census +
## mod.tree_adult6: A.con_curt + A.het_curt + (1 | quad.unique) + (1 | sp.)
## mod.tree_adult4: cbind(survs, deaths) ~ site + pesticide * exclosure * census +
## mod.tree_adult4: A.con_curt + A.het_curt + (1 | quad.unique) + (1 | sp.)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## mod.tree_adult5 17 5272.3 5371.9 -2619.1 5238.3
## mod.tree_adult6 26 5254.2 5406.6 -2601.1 5202.2 36.009 9 3.951e-05
## mod.tree_adult4 38 5268.3 5491.0 -2596.2 5192.3 9.938 12 0.6214
##
## mod.tree_adult5
## mod.tree_adult6 ***
## mod.tree_adult4
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# without neighoring trees
mod.tree_adult7 <- glmer(cbind(survs, deaths) ~ site + pesticide + exclosure + census +
(1|quad.unique) + (1|sp.),
data=subset(pestdat.ag_adult, growth.form=='tree'), family=binomial,
control=glmerControl(optimizer="optimx",
optCtrl=list(method=c('bobyqa', 'Nelder-Mead'))))
# nested random effects
mod.tree_adult8 <- glmer(cbind(survs, deaths) ~ pesticide + exclosure + census +
A.con_curt + A.het_curt +
(1|plot/quadrat) + (1|sp.),
data=subset(pestdat.ag_adult, growth.form=='tree'), family=binomial,
control=glmerControl(optimizer="optimx",
optCtrl=list(method=c('bobyqa', 'Nelder-Mead'))))
anova(mod.tree_adult1, mod.tree_adult8)## Data: subset(pestdat.ag_adult, growth.form == "tree")
## Models:
## mod.tree_adult8: cbind(survs, deaths) ~ pesticide + exclosure + census + A.con_curt +
## mod.tree_adult8: A.het_curt + (1 | plot/quadrat) + (1 | sp.)
## mod.tree_adult1: cbind(survs, deaths) ~ site + pesticide + exclosure + census +
## mod.tree_adult1: A.con_curt + A.het_curt + (1 | quad.unique) + (1 | sp.)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## mod.tree_adult8 13 5269.0 5345.2 -2621.5 5243.0
## mod.tree_adult1 14 5269.1 5351.2 -2620.6 5241.1 1.8897 1 0.1692mod.tree_adult9 <- glmer(cbind(survs, deaths) ~ pesticide + exclosure +
A.con_curt + A.het_curt +
(1|plot/quadrat) + (1|plot/census) + (1|sp.),
data=subset(pestdat.ag_adult, growth.form=='tree'), family=binomial,
control=glmerControl(optimizer="optimx",
optCtrl=list(method=c('bobyqa', 'Nelder-Mead'))))## Warning in optimx.check(par, optcfg$ufn, optcfg$ugr, optcfg$uhess, lower, : Parameters or bounds appear to have different scalings.
## This can cause poor performance in optimization.
## It is important for derivative free methods like BOBYQA, UOBYQA, NEWUOA.mod.tree_adult10 <- glmer(cbind(survs, deaths) ~ pesticide + exclosure + census +
A.con_curt + A.het_curt +
(1|plot/pesticide/quadrat) + (1|plot/census) + (1|sp.),
data=subset(pestdat.ag_adult, growth.form=='tree'), family=binomial,
control=glmerControl(optimizer="optimx",
optCtrl=list(method=c('bobyqa', 'Nelder-Mead'))))## Warning in optimx.check(par, optcfg$ufn, optcfg$ugr, optcfg$uhess, lower, : Parameters or bounds appear to have different scalings.
## This can cause poor performance in optimization.
## It is important for derivative free methods like BOBYQA, UOBYQA, NEWUOA.## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, : Model is nearly unidentifiable: large eigenvalue ratio
## - Rescale variables?mod.tree_adult11 <- glmer(cbind(survs, deaths) ~ pesticide + exclosure + census +
A.con_curt + A.het_curt +
(1|quad.unique) + (1|sp.),
data=subset(pestdat.ag_adult, growth.form=='tree'), family=binomial,
control=glmerControl(optimizer="optimx",
optCtrl=list(method=c('bobyqa', 'Nelder-Mead'))))
anova(mod.tree_adult1, mod.tree_adult8, mod.tree_adult11)## Data: subset(pestdat.ag_adult, growth.form == "tree")
## Models:
## mod.tree_adult11: cbind(survs, deaths) ~ pesticide + exclosure + census + A.con_curt +
## mod.tree_adult11: A.het_curt + (1 | quad.unique) + (1 | sp.)
## mod.tree_adult8: cbind(survs, deaths) ~ pesticide + exclosure + census + A.con_curt +
## mod.tree_adult8: A.het_curt + (1 | plot/quadrat) + (1 | sp.)
## mod.tree_adult1: cbind(survs, deaths) ~ site + pesticide + exclosure + census +
## mod.tree_adult1: A.con_curt + A.het_curt + (1 | quad.unique) + (1 | sp.)
## Df AIC BIC logLik deviance Chisq Chi Df
## mod.tree_adult11 12 5276.4 5346.8 -2626.2 5252.4
## mod.tree_adult8 13 5269.0 5345.2 -2621.5 5243.0 9.4309 1
## mod.tree_adult1 14 5269.1 5351.2 -2620.6 5241.1 1.8897 1
## Pr(>Chisq)
## mod.tree_adult11
## mod.tree_adult8 0.002134 **
## mod.tree_adult1 0.169234
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1mod.tree_adult12 <- glmer(cbind(survs, deaths) ~ site + pesticide*exclosure*census +
A.con_curt + A.het_curt +
(1|quad.unique) + (1|sp.),
data=subset(pestdat.ag_adult, growth.form=='tree'), family=binomial,
control=glmerControl(optimizer="optimx",
optCtrl=list(method=c('bobyqa', 'Nelder-Mead'))))
mod.tree_adult13 <- glmer(cbind(survs, deaths) ~ site + (pesticide+exclosure)*census +
A.con_curt + A.het_curt +
(1|quad.unique) + (1|sp.),
data=subset(pestdat.ag_adult, growth.form=='tree'), family=binomial,
control=glmerControl(optimizer="optimx",
optCtrl=list(method=c('bobyqa', 'Nelder-Mead'))))
mod.tree_adult14 <- glmer(cbind(survs, deaths) ~ site + pesticide*census + exclosure+
A.con_curt + A.het_curt +
(1|quad.unique) + (1|sp.),
data=subset(pestdat.ag_adult, growth.form=='tree'), family=binomial,
control=glmerControl(optimizer="optimx",
optCtrl=list(method=c('bobyqa', 'Nelder-Mead'))))
mod.tree_adult15 <- glmer(cbind(survs, deaths) ~ pesticide*census + exclosure+
A.con_curt + A.het_curt +
(1|plot/quadrat) + (1|sp.),
data=subset(pestdat.ag_adult, growth.form=='tree'), family=binomial,
control=glmerControl(optimizer="optimx",
optCtrl=list(method=c('bobyqa', 'Nelder-Mead'))))
anova(mod.tree_adult1, mod.tree_adult8, mod.tree_adult13, mod.tree_adult14, mod.tree_adult15)## Data: subset(pestdat.ag_adult, growth.form == "tree")
## Models:
## mod.tree_adult8: cbind(survs, deaths) ~ pesticide + exclosure + census + A.con_curt +
## mod.tree_adult8: A.het_curt + (1 | plot/quadrat) + (1 | sp.)
## mod.tree_adult1: cbind(survs, deaths) ~ site + pesticide + exclosure + census +
## mod.tree_adult1: A.con_curt + A.het_curt + (1 | quad.unique) + (1 | sp.)
## mod.tree_adult15: cbind(survs, deaths) ~ pesticide * census + exclosure + A.con_curt +
## mod.tree_adult15: A.het_curt + (1 | plot/quadrat) + (1 | sp.)
## mod.tree_adult14: cbind(survs, deaths) ~ site + pesticide * census + exclosure +
## mod.tree_adult14: A.con_curt + A.het_curt + (1 | quad.unique) + (1 | sp.)
## mod.tree_adult13: cbind(survs, deaths) ~ site + (pesticide + exclosure) * census +
## mod.tree_adult13: A.con_curt + A.het_curt + (1 | quad.unique) + (1 | sp.)
## Df AIC BIC logLik deviance Chisq Chi Df
## mod.tree_adult8 13 5269.0 5345.2 -2621.5 5243.0
## mod.tree_adult1 14 5269.1 5351.2 -2620.6 5241.1 1.8897 1
## mod.tree_adult15 22 5251.3 5380.2 -2603.7 5207.3 33.7936 8
## mod.tree_adult14 23 5251.5 5386.3 -2602.8 5205.5 1.8032 1
## mod.tree_adult13 26 5254.2 5406.6 -2601.1 5202.2 3.2734 3
## Pr(>Chisq)
## mod.tree_adult8
## mod.tree_adult1 0.1692
## mod.tree_adult15 4.428e-05 ***
## mod.tree_adult14 0.1793
## mod.tree_adult13 0.3514
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1summary(mod.tree_adult14, correlation=FALSE)## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: cbind(survs, deaths) ~ site + pesticide * census + exclosure +
## A.con_curt + A.het_curt + (1 | quad.unique) + (1 | sp.)
## Data: subset(pestdat.ag_adult, growth.form == "tree")
## Control:
## glmerControl(optimizer = "optimx", optCtrl = list(method = c("bobyqa",
## "Nelder-Mead")))
##
## AIC BIC logLik deviance df.resid
## 5251.5 5386.3 -2602.8 5205.5 2568
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -4.6563 -0.5954 0.3777 0.6731 2.3912
##
## Random effects:
## Groups Name Variance Std.Dev.
## quad.unique (Intercept) 0.2069 0.4548
## sp. (Intercept) 0.5349 0.7314
## Number of obs: 2591, groups: quad.unique, 239; sp., 20
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.53085 0.70253 -0.756 0.449879
## site1 0.01310 0.06051 0.216 0.828655
## site2 -0.17776 0.05713 -3.111 0.001862 **
## pesticideF 0.73671 0.15200 4.847 1.25e-06 ***
## pesticideI 0.11356 0.15168 0.749 0.454041
## pesticideC 0.03071 0.15330 0.200 0.841225
## census16sp -1.08980 0.22560 -4.831 1.36e-06 ***
## census17fa 0.34061 0.13116 2.597 0.009405 **
## census17sp -0.46765 0.12700 -3.682 0.000231 ***
## exclosure1 0.09676 0.08302 1.166 0.243816
## A.con_curt -0.03045 0.01342 -2.269 0.023293 *
## A.het_curt 0.13529 0.03718 3.639 0.000274 ***
## pesticideF:census16sp -0.89672 0.31745 -2.825 0.004731 **
## pesticideI:census16sp -0.51417 0.33445 -1.537 0.124198
## pesticideC:census16sp -0.04847 0.33054 -0.147 0.883426
## pesticideF:census17fa -0.92802 0.17629 -5.264 1.41e-07 ***
## pesticideI:census17fa -0.24958 0.18305 -1.363 0.172740
## pesticideC:census17fa -0.34055 0.18703 -1.821 0.068632 .
## pesticideF:census17sp -0.61059 0.17587 -3.472 0.000517 ***
## pesticideI:census17sp -0.20690 0.18029 -1.148 0.251123
## pesticideC:census17sp -0.22216 0.18457 -1.204 0.228699
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1##
## Correlation matrix not shown by default, as p = 21 > 12.
## Use print(x, correlation=TRUE) or
## vcov(x) if you need it## select the approprate model and test the sensitivity
mod.tree_adult.diags <- augment(mod.tree_adult1)
qplot(data=mod.tree_adult.diags, x=.fitted, y=.resid, geom=c('point', 'smooth')) +
geom_hline(yintercept=0, linetype='dotted')## `geom_smooth()` using method = 'gam'
qplot(data = mod.tree_adult.diags, x= pesticide, y=.wtres, geom="boxplot", colour=exclosure) 
library(car)##
## Attaching package: 'car'## The following object is masked from 'package:dplyr':
##
## recodeqqPlot(resid(mod.tree_adult1))
qqPlot(resid(mod.tree_adult14))
library(sjPlot)## Visit http://strengejacke.de/sjPlot for package-vignettes.##
## Attaching package: 'sjPlot'## The following objects are masked from 'package:cowplot':
##
## plot_grid, save_plotsjp.lmer(mod.tree_adult1, type='fe')
sjp.lmer(mod.tree_adult1, type='re', sort.est="sort.all" )## Sorting each group of random intercept ('sort.all') is not possible when 'facet.grid = TRUE'.## Plotting random effects...## Sorting each group of random intercept ('sort.all') is not possible when 'facet.grid = TRUE'.## Plotting random effects...
summary(mod.tree_adult1, correlation=FALSE)## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: cbind(survs, deaths) ~ site + pesticide + exclosure + census +
## A.con_curt + A.het_curt + (1 | quad.unique) + (1 | sp.)
## Data: subset(pestdat.ag_adult, growth.form == "tree")
## Control:
## glmerControl(optimizer = "optimx", optCtrl = list(method = c("bobyqa",
## "Nelder-Mead")))
##
## AIC BIC logLik deviance df.resid
## 5269.1 5351.2 -2620.6 5241.1 2577
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -4.6535 -0.5888 0.3808 0.6746 2.3480
##
## Random effects:
## Groups Name Variance Std.Dev.
## quad.unique (Intercept) 0.2088 0.4570
## sp. (Intercept) 0.5085 0.7131
## Number of obs: 2591, groups: quad.unique, 239; sp., 20
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.31548 0.69975 -0.451 0.652105
## site1 0.01390 0.06058 0.229 0.818543
## site2 -0.17552 0.05721 -3.068 0.002155 **
## pesticideF 0.19337 0.11067 1.747 0.080599 .
## pesticideI -0.04799 0.11446 -0.419 0.675018
## pesticideC -0.13903 0.11502 -1.209 0.226767
## exclosure1 0.09986 0.08308 1.202 0.229380
## census16sp -1.45024 0.12070 -12.015 < 2e-16 ***
## census17fa -0.06942 0.06558 -1.058 0.289829
## census17sp -0.73059 0.06433 -11.356 < 2e-16 ***
## A.con_curt -0.03002 0.01341 -2.239 0.025128 *
## A.het_curt 0.13482 0.03723 3.621 0.000293 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1## shrub
mod.shrub_adult1 <- glmer(cbind(survs, deaths) ~ site + pesticide + exclosure + census + A.sum_curt +
(1|quad.unique) + (1|sp.),
data=subset(pestdat.ag_adult, growth.form=='shrub'), family=binomial,
control=glmerControl(optimizer="optimx",
optCtrl=list(method=c('bobyqa', 'Nelder-Mead'))))## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control
## $checkConv, : unable to evaluate scaled gradient## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control
## $checkConv, : Model failed to converge: degenerate Hessian with 1 negative
## eigenvalues## model without sum area of adult neigboring trees
mod.shrub_adult2 <- glmer(cbind(survs, deaths) ~ pesticide + exclosure + census + (1|quad.unique) + (1|sp.),
data=subset(pestdat.ag_adult, growth.form=='shrub'), family=binomial,
control=glmerControl(optimizer="optimx",
optCtrl=list(method=c('bobyqa', 'Nelder-Mead'))))
mod.shrub_adult3 <- glmer(cbind(survs, deaths) ~ pesticide + exclosure + census + (1|plot/quadrat) + (1|sp.),
data=subset(pestdat.ag_adult, growth.form=='shrub'), family=binomial,
control=glmerControl(optimizer="optimx",
optCtrl=list(method=c('bobyqa', 'Nelder-Mead'))))
mod.shrub_adult4 <- glmer(cbind(survs, deaths) ~ pesticide + exclosure + census + A.sum_curt +
(1|plot/quadrat) + (1|sp.),
data=subset(pestdat.ag_adult, growth.form=='shrub'), family=binomial,
control=glmerControl(optimizer="optimx",
optCtrl=list(method=c('bobyqa', 'Nelder-Mead'))))## Warning in optimx.check(par, optcfg$ufn, optcfg$ugr, optcfg$uhess, lower, : Parameters or bounds appear to have different scalings.
## This can cause poor performance in optimization.
## It is important for derivative free methods like BOBYQA, UOBYQA, NEWUOA.mod.shrub_adult5 <- glmer(cbind(survs, deaths) ~ site + pesticide + exclosure + census +
(1|quad.unique) + (1|sp.),
data=subset(pestdat.ag_adult, growth.form=='shrub'), family=binomial,
control=glmerControl(optimizer="optimx",
optCtrl=list(method=c('bobyqa', 'Nelder-Mead'))))
mod.shrub_adult6 <- glmer(cbind(survs, deaths) ~ site + pesticide*census + exclosure +
(1|quad.unique) + (1|sp.),
data=subset(pestdat.ag_adult, growth.form=='shrub'), family=binomial,
control=glmerControl(optimizer="optimx",
optCtrl=list(method=c('bobyqa', 'Nelder-Mead'))))## Warning in optimx.check(par, optcfg$ufn, optcfg$ugr, optcfg$uhess, lower, : Parameters or bounds appear to have different scalings.
## This can cause poor performance in optimization.
## It is important for derivative free methods like BOBYQA, UOBYQA, NEWUOA.## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, : Model is nearly unidentifiable: large eigenvalue ratio
## - Rescale variables?anova(mod.shrub_adult1, mod.shrub_adult3, mod.shrub_adult4, mod.shrub_adult5, mod.shrub_adult6)## Data: subset(pestdat.ag_adult, growth.form == "shrub")
## Models:
## mod.shrub_adult3: cbind(survs, deaths) ~ pesticide + exclosure + census + (1 |
## mod.shrub_adult3: plot/quadrat) + (1 | sp.)
## mod.shrub_adult4: cbind(survs, deaths) ~ pesticide + exclosure + census + A.sum_curt +
## mod.shrub_adult4: (1 | plot/quadrat) + (1 | sp.)
## mod.shrub_adult5: cbind(survs, deaths) ~ site + pesticide + exclosure + census +
## mod.shrub_adult5: (1 | quad.unique) + (1 | sp.)
## mod.shrub_adult1: cbind(survs, deaths) ~ site + pesticide + exclosure + census +
## mod.shrub_adult1: A.sum_curt + (1 | quad.unique) + (1 | sp.)
## mod.shrub_adult6: cbind(survs, deaths) ~ site + pesticide * census + exclosure +
## mod.shrub_adult6: (1 | quad.unique) + (1 | sp.)
## Df AIC BIC logLik deviance Chisq Chi Df
## mod.shrub_adult3 11 486.90 546.98 -232.45 464.90
## mod.shrub_adult4 12 488.18 553.73 -232.09 464.18 0.7110 1
## mod.shrub_adult5 12 486.44 551.98 -231.22 462.44 1.7487 0
## mod.shrub_adult1 13 487.26 558.27 -230.63 461.26 1.1793 1
## mod.shrub_adult6 21 496.75 611.46 -227.38 454.75 6.5059 8
## Pr(>Chisq)
## mod.shrub_adult3
## mod.shrub_adult4 0.3991
## mod.shrub_adult5 <2e-16 ***
## mod.shrub_adult1 0.2775
## mod.shrub_adult6 0.5908
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# add sum of neighoring trees make no sense
# model sensitivity tests
mod.shrub_adult.diags <- augment(mod.shrub_adult5)
qplot(data=mod.shrub_adult.diags, x=.fitted, y=.resid, geom=c('point', 'smooth')) +
geom_hline(yintercept=0, linetype='dotted')## `geom_smooth()` using method = 'gam'
qplot(data = mod.shrub_adult.diags, x= pesticide, y=.wtres, geom="boxplot", colour=exclosure) 
qqPlot(resid(mod.shrub_adult5))
sjp.lmer(mod.shrub_adult5, type='fe')
sjp.lmer(mod.shrub_adult5, type='re', sort.est="sort.all" )## Sorting each group of random intercept ('sort.all') is not possible when 'facet.grid = TRUE'.
## Plotting random effects...## Sorting each group of random intercept ('sort.all') is not possible when 'facet.grid = TRUE'.## Plotting random effects...
summary(mod.shrub_adult5, correlation=FALSE)## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: cbind(survs, deaths) ~ site + pesticide + exclosure + census +
## (1 | quad.unique) + (1 | sp.)
## Data: subset(pestdat.ag_adult, growth.form == "shrub")
## Control:
## glmerControl(optimizer = "optimx", optCtrl = list(method = c("bobyqa",
## "Nelder-Mead")))
##
## AIC BIC logLik deviance df.resid
## 486.4 552.0 -231.2 462.4 1729
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -12.2932 0.0906 0.1282 0.1877 0.5913
##
## Random effects:
## Groups Name Variance Std.Dev.
## quad.unique (Intercept) 0.35332 0.5944
## sp. (Intercept) 0.09279 0.3046
## Number of obs: 1741, groups: quad.unique, 195; sp., 19
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 4.86943 0.56128 8.676 < 2e-16 ***
## site1 -0.15280 0.20412 -0.749 0.454108
## site2 0.33838 0.23371 1.448 0.147657
## pesticideF -0.69626 0.40450 -1.721 0.085197 .
## pesticideI -0.17234 0.44019 -0.392 0.695412
## pesticideC 0.02797 0.45794 0.061 0.951289
## exclosure1 0.56311 0.31036 1.814 0.069616 .
## census16sp -1.65148 0.45466 -3.632 0.000281 ***
## census17fa 0.06338 0.58204 0.109 0.913291
## census17sp -1.11219 0.47702 -2.332 0.019724 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Robi, I have one question here. For the tree survival analysis, the model considering interation between pesticide and census fits better than those without the interaction. Should we keep the interacton term? (See code: anova(mod.tree_adult1, mod.tree_adult8, mod.tree_adult13, mod.tree_adult14, mod.tree_adult15) Line 298) If so, how to interpret the results? I mean, several interactions between fungicide and censuses are negative. Is that mean they compared with the water treatment in first census?
Density dependence
##------------------------- Density dependence-------------------------##
## For all tree species together
ggplot(data=subset(pestdat.ag_adult, growth.form=='tree'),
aes(x=total, y=survs/total, colour=pesticide)) +
geom_jitter()+
geom_smooth(aes(weight = total), method='glm', method.args=list(family="binomial")) +
facet_wrap(~exclosure)
# hard to transform the density of seendlings to normal distributed
hist(pestdat.ag_adult$total)
hist(log10(pestdat.ag_adult$total))
pestdat.ag_adult$dens <- (pestdat.ag_adult$total - mean(pestdat.ag_adult$total))/(2*sd(pestdat.ag_adult$total))
mod.tree.ddst_adult1 <- glmer(cbind(survs, deaths) ~ site + pesticide + exclosure + census + A.con_curt + A.het_curt + dens + (1|quad.unique) + (1|sp.),
data=subset(pestdat.ag_adult, growth.form=='tree'), family=binomial,
control=glmerControl(optimizer="optimx",
optCtrl=list(method=c('bobyqa', 'Nelder-Mead'))))
mod.tree.ddst_adult2 <- glmer(cbind(survs, deaths) ~ pesticide + exclosure + census + A.con_curt + A.het_curt + dens + (1|plot/quadrat) + (1|sp.),
data=subset(pestdat.ag_adult, growth.form=='tree'), family=binomial,
control=glmerControl(optimizer="optimx",
optCtrl=list(method=c('bobyqa', 'Nelder-Mead'))))
anova(mod.tree_adult8, mod.tree.ddst_adult1, mod.tree.ddst_adult2)## Data: subset(pestdat.ag_adult, growth.form == "tree")
## Models:
## mod.tree_adult8: cbind(survs, deaths) ~ pesticide + exclosure + census + A.con_curt +
## mod.tree_adult8: A.het_curt + (1 | plot/quadrat) + (1 | sp.)
## mod.tree.ddst_adult2: cbind(survs, deaths) ~ pesticide + exclosure + census + A.con_curt +
## mod.tree.ddst_adult2: A.het_curt + dens + (1 | plot/quadrat) + (1 | sp.)
## mod.tree.ddst_adult1: cbind(survs, deaths) ~ site + pesticide + exclosure + census +
## mod.tree.ddst_adult1: A.con_curt + A.het_curt + dens + (1 | quad.unique) + (1 |
## mod.tree.ddst_adult1: sp.)
## Df AIC BIC logLik deviance Chisq Chi Df
## mod.tree_adult8 13 5269.0 5345.2 -2621.5 5243.0
## mod.tree.ddst_adult2 14 5265.7 5347.7 -2618.8 5237.7 5.3047 1
## mod.tree.ddst_adult1 15 5266.2 5354.1 -2618.1 5236.2 1.4530 1
## Pr(>Chisq)
## mod.tree_adult8
## mod.tree.ddst_adult2 0.02127 *
## mod.tree.ddst_adult1 0.22805
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1mod.tree.ddst_adult.diags <- augment(mod.tree.ddst_adult1)
qplot(data=mod.tree.ddst_adult.diags, x=.fitted, y=.resid, geom=c('point', 'smooth')) +
geom_hline(yintercept=0, linetype='dotted')## `geom_smooth()` using method = 'gam'
qplot(data =mod.tree.ddst_adult.diags, x= pesticide, y=.wtres, geom="boxplot", colour=exclosure) 
qqPlot(resid(mod.tree.ddst_adult1))
sjp.lmer(mod.tree.ddst_adult1, type='fe')
sjp.lmer(mod.tree.ddst_adult1, type='re', sort.est="sort.all" )## Sorting each group of random intercept ('sort.all') is not possible when 'facet.grid = TRUE'.## Plotting random effects...## Sorting each group of random intercept ('sort.all') is not possible when 'facet.grid = TRUE'.## Plotting random effects...
summary(mod.tree.ddst_adult1, correlation=FALSE)## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: cbind(survs, deaths) ~ site + pesticide + exclosure + census +
## A.con_curt + A.het_curt + dens + (1 | quad.unique) + (1 | sp.)
## Data: subset(pestdat.ag_adult, growth.form == "tree")
## Control:
## glmerControl(optimizer = "optimx", optCtrl = list(method = c("bobyqa",
## "Nelder-Mead")))
##
## AIC BIC logLik deviance df.resid
## 5266.2 5354.1 -2618.1 5236.2 2576
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -4.2386 -0.6019 0.3719 0.6590 2.3434
##
## Random effects:
## Groups Name Variance Std.Dev.
## quad.unique (Intercept) 0.2089 0.4570
## sp. (Intercept) 0.4982 0.7058
## Number of obs: 2591, groups: quad.unique, 239; sp., 20
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.370323 0.700573 -0.529 0.59708
## site1 0.009565 0.060692 0.158 0.87477
## site2 -0.155523 0.057960 -2.683 0.00729 **
## pesticideF 0.212304 0.111082 1.911 0.05597 .
## pesticideI -0.049563 0.114523 -0.433 0.66517
## pesticideC -0.148329 0.115166 -1.288 0.19776
## exclosure1 0.104011 0.083177 1.250 0.21113
## census16sp -1.516264 0.124482 -12.181 < 2e-16 ***
## census17fa -0.083965 0.066092 -1.270 0.20393
## census17sp -0.757571 0.065556 -11.556 < 2e-16 ***
## A.con_curt -0.027493 0.013494 -2.037 0.04161 *
## A.het_curt 0.138247 0.037296 3.707 0.00021 ***
## dens -0.091498 0.041339 -2.213 0.02687 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1##
## Correlation matrix not shown by default, as p = 13 > 12.
## Use print(x, correlation=TRUE) or
## vcov(x) if you need it##------------------------- individual species -------------------------------##
## And now pulling the species apart.
## First need to filter down to species that can be analysed - need enough
## reps (lets say 10) and variation in abundance (at least > 5
sp.counts <- summarise(group_by(pestdat.ag_adult, sp.),
n.quadrat=length(total), max.dens=max(total), min.dens=min(total))
sp.counts$range <- sp.counts$max.dens - sp.counts$min.dens
sp.sel <- sp.counts$sp.[sp.counts$n.quadrat > 10 & sp.counts$range > 5]
sp.counts[sp.counts$sp. %in% sp.sel,]## # A tibble: 8 × 5
## sp. n.quadrat max.dens min.dens range
## <fctr> <int> <int> <int> <int>
## 1 ABNE 46 12 1 11
## 2 ACBA 260 19 1 18
## 3 ACPS 296 12 1 11
## 4 ACTE 147 8 1 7
## 5 FRMA 765 23 1 22
## 6 PIKO 75 8 1 7
## 7 QUMO 52 8 1 7
## 8 TIAM 673 22 1 21ggplot(data=subset(pestdat.ag_adult, sp. %in% sp.sel),
aes(x=total, y=survs/total, colour=pesticide)) +
geom_jitter()+
geom_smooth(aes(weight = total), method='glm',
method.args=list(family="binomial"), fullrange=T) +
facet_wrap(~ sp., scales='free')
# it seems only eight species have enough data, but conference intervals of ACMO and QUMO are too huge
## quick plot of only the species for which there appear to be enough data
# sp.sel2 <- c('ABNE', 'ACBA', 'ACMO', 'ACPS', 'ACTE', 'FRMA', 'PIKO', 'QUMO', 'TIAM')
sp.sel2 <- c('ACBA', 'ACPS', 'ACTE', 'FRMA', 'PIKO', 'TIAM')
ggplot(data=subset(pestdat.ag_adult, sp. %in% sp.sel2),
aes(x=total, y=survs/total, colour=pesticide)) +
geom_jitter()+
geom_smooth(aes(weight = total), method='glm',
method.args=list(family="binomial"), fullrange=T) +
facet_wrap(~ sp. + exclosure, scales='free')
nddmods <- sapply(sp.sel2, function(sp){
print(sp)
spdat <- filter(pestdat.ag_adult, sp. == sp)
## to increase ability to fit models, standardise total
spdat$dens <- (spdat$total - mean(spdat$total))/(2*sd(spdat$total))
spdat <- droplevels(spdat)
mod <- glmer(cbind(survs, deaths) ~ site + pesticide + exclosure + census +
A.con_curt + A.het_curt + dens + (1|quad.unique),
data=spdat, family=binomial,
control=glmerControl(optimizer="optimx",
optCtrl=list(method=c('bobyqa', 'Nelder-Mead')))
)
}, simplify=FALSE)## [1] "ACBA"
## [1] "ACPS"
## [1] "ACTE"## Warning in optimx.check(par, optcfg$ufn, optcfg$ugr, optcfg$uhess, lower, : Parameters or bounds appear to have different scalings.
## This can cause poor performance in optimization.
## It is important for derivative free methods like BOBYQA, UOBYQA, NEWUOA.## [1] "FRMA"
## [1] "PIKO"
## [1] "TIAM"## diagnostic plots
indmods.resids <- lapply(nddmods, augment)
indmods.resids <- do.call('rbind', mapply(function(spnm, dat) {
dat$species <- spnm
return(dat)
}, spnm=as.list(names(indmods.resids)), dat=indmods.resids,
SIMPLIFY=FALSE))
ggplot(data=indmods.resids, aes(x=.fitted, y=.wtres, colour=pesticide)) +
geom_jitter() + geom_smooth() + facet_wrap(~species, scales='free') +
geom_hline(yintercept=0, linetype='dashed')## `geom_smooth()` using method = 'loess'
library(lattice)
lapply(nddmods, function(mod){
dotplot(ranef(mod, condVar=TRUE))})## $ACBA
## $ACBA$quad.unique
##
##
## $ACPS
## $ACPS$quad.unique
##
##
## $ACTE
## $ACTE$quad.unique
##
##
## $FRMA
## $FRMA$quad.unique
##
##
## $PIKO
## $PIKO$quad.unique
##
##
## $TIAM
## $TIAM$quad.unique
lapply(nddmods, function(x) summary(x, correlation=FALSE))## $ACBA
## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: cbind(survs, deaths) ~ site + pesticide + exclosure + census +
## A.con_curt + A.het_curt + dens + (1 | quad.unique)
## Data: spdat
## Control:
## glmerControl(optimizer = "optimx", optCtrl = list(method = c("bobyqa",
## "Nelder-Mead")))
##
## AIC BIC logLik deviance df.resid
## 383.7 433.5 -177.8 355.7 246
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.9506 -0.6146 0.3482 0.5208 1.3918
##
## Random effects:
## Groups Name Variance Std.Dev.
## quad.unique (Intercept) 1.378 1.174
## Number of obs: 260, groups: quad.unique, 114
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 3.43543 4.08147 0.842 0.39995
## site2 -1.36234 0.51243 -2.659 0.00785 **
## site3 -1.04899 0.61074 -1.718 0.08588 .
## pesticideF -0.27307 0.52806 -0.517 0.60508
## pesticideI -0.38875 0.52778 -0.737 0.46138
## pesticideC -0.05273 0.56505 -0.093 0.92565
## exclosure1 -0.59540 0.43482 -1.369 0.17091
## census16sp -0.70149 0.60670 -1.156 0.24759
## census17fa 0.42763 0.36927 1.158 0.24685
## census17sp -1.01130 0.35893 -2.817 0.00484 **
## A.con_curt -0.01353 0.32369 -0.042 0.96667
## A.het_curt -0.04527 0.22333 -0.203 0.83937
## dens -0.44753 0.28677 -1.561 0.11861
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1##
## Correlation matrix not shown by default, as p = 13 > 12.
## Use print(...., correlation=TRUE) or
## vcov(....) if you need it##
## $ACPS
## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: cbind(survs, deaths) ~ site + pesticide + exclosure + census +
## A.con_curt + A.het_curt + dens + (1 | quad.unique)
## Data: spdat
## Control:
## glmerControl(optimizer = "optimx", optCtrl = list(method = c("bobyqa",
## "Nelder-Mead")))
##
## AIC BIC logLik deviance df.resid
## 370.1 421.7 -171.0 342.1 282
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -8.3717 0.0843 0.2185 0.4799 1.1045
##
## Random effects:
## Groups Name Variance Std.Dev.
## quad.unique (Intercept) 0.4266 0.6531
## Number of obs: 296, groups: quad.unique, 159
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 4.60321 2.84794 1.616 0.10602
## site2 -0.08561 0.31637 -0.271 0.78671
## site3 0.31949 0.41168 0.776 0.43772
## pesticideF 0.08305 0.35728 0.232 0.81618
## pesticideI 0.12886 0.38874 0.332 0.74029
## pesticideC -0.08517 0.39967 -0.213 0.83125
## exclosure1 -0.22560 0.33685 -0.670 0.50302
## census16sp -2.49543 1.27919 -1.951 0.05108 .
## census17fa -3.20721 1.02650 -3.124 0.00178 **
## census17sp -1.67945 1.10957 -1.514 0.13012
## A.con_curt -0.11535 0.15830 -0.729 0.46617
## A.het_curt 0.04720 0.15061 0.313 0.75397
## dens 0.22956 0.22045 1.041 0.29772
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1##
## Correlation matrix not shown by default, as p = 13 > 12.
## Use print(...., correlation=TRUE) or
## vcov(....) if you need it##
## $ACTE
## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: cbind(survs, deaths) ~ site + pesticide + exclosure + census +
## A.con_curt + A.het_curt + dens + (1 | quad.unique)
## Data: spdat
## Control:
## glmerControl(optimizer = "optimx", optCtrl = list(method = c("bobyqa",
## "Nelder-Mead")))
##
## AIC BIC logLik deviance df.resid
## 188.1 227.0 -81.1 162.1 134
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -4.8905 -0.1223 0.3494 0.5458 1.2432
##
## Random effects:
## Groups Name Variance Std.Dev.
## quad.unique (Intercept) 0 0
## Number of obs: 147, groups: quad.unique, 59
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -5.62217 4.29189 -1.310 0.1902
## site2 -0.15382 0.43979 -0.350 0.7265
## pesticideF -0.73375 0.55433 -1.324 0.1856
## pesticideI -0.63145 0.52160 -1.211 0.2260
## pesticideC -0.51642 0.56344 -0.917 0.3594
## exclosure1 -0.11001 0.55640 -0.198 0.8433
## census16sp -1.48944 0.62940 -2.366 0.0180 *
## census17fa 0.44397 0.45710 0.971 0.3314
## census17sp 0.20571 0.47039 0.437 0.6619
## A.con_curt -0.47809 0.28447 -1.681 0.0928 .
## A.het_curt 0.54951 0.23876 2.301 0.0214 *
## dens -0.05784 0.27089 -0.214 0.8309
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## convergence code: 0
## Parameters or bounds appear to have different scalings.
## This can cause poor performance in optimization.
## It is important for derivative free methods like BOBYQA, UOBYQA, NEWUOA.
##
##
## $FRMA
## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: cbind(survs, deaths) ~ site + pesticide + exclosure + census +
## A.con_curt + A.het_curt + dens + (1 | quad.unique)
## Data: spdat
## Control:
## glmerControl(optimizer = "optimx", optCtrl = list(method = c("bobyqa",
## "Nelder-Mead")))
##
## AIC BIC logLik deviance df.resid
## 1680.9 1745.8 -826.4 1652.9 751
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.9632 -0.6714 0.3874 0.6946 1.9496
##
## Random effects:
## Groups Name Variance Std.Dev.
## quad.unique (Intercept) 0.2132 0.4617
## Number of obs: 765, groups: quad.unique, 230
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.050470 0.929449 0.054 0.9567
## site2 -0.354717 0.175897 -2.017 0.0437 *
## site3 -0.180762 0.179258 -1.008 0.3133
## pesticideF 0.062224 0.163599 0.380 0.7037
## pesticideI -0.301626 0.162044 -1.861 0.0627 .
## pesticideC -0.139590 0.165317 -0.844 0.3985
## exclosure1 0.251150 0.120128 2.091 0.0366 *
## census16sp -1.560813 0.176617 -8.837 < 2e-16 ***
## census17fa 0.310369 0.136724 2.270 0.0232 *
## census17sp -0.495096 0.108348 -4.570 4.89e-06 ***
## A.con_curt 0.006949 0.016989 0.409 0.6825
## A.het_curt 0.114272 0.052438 2.179 0.0293 *
## dens -0.179276 0.120689 -1.485 0.1374
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1##
## Correlation matrix not shown by default, as p = 13 > 12.
## Use print(...., correlation=TRUE) or
## vcov(....) if you need it##
## $PIKO
## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: cbind(survs, deaths) ~ site + pesticide + exclosure + census +
## A.con_curt + A.het_curt + dens + (1 | quad.unique)
## Data: spdat
## Control:
## glmerControl(optimizer = "optimx", optCtrl = list(method = c("bobyqa",
## "Nelder-Mead")))
##
## AIC BIC logLik deviance df.resid
## 146.7 179.2 -59.4 118.7 61
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -4.8041 -0.5187 0.3436 0.6499 1.1901
##
## Random effects:
## Groups Name Variance Std.Dev.
## quad.unique (Intercept) 0.6002 0.7748
## Number of obs: 75, groups: quad.unique, 43
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -6.89039 6.98760 -0.986 0.32409
## site2 -0.17103 0.83838 -0.204 0.83836
## site3 -0.80554 0.95573 -0.843 0.39931
## pesticideF 2.61863 0.95197 2.751 0.00595 **
## pesticideI 0.66617 0.85409 0.780 0.43541
## pesticideC 0.76070 0.99930 0.761 0.44651
## exclosure1 0.24621 0.62163 0.396 0.69206
## census16sp -1.45947 1.96659 -0.742 0.45801
## census17fa 0.53421 0.71735 0.745 0.45645
## census17sp -0.46806 0.59740 -0.784 0.43334
## A.con_curt 0.62626 0.42384 1.478 0.13952
## A.het_curt 0.05408 0.24498 0.221 0.82529
## dens 0.07058 0.51695 0.136 0.89140
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1##
## Correlation matrix not shown by default, as p = 13 > 12.
## Use print(...., correlation=TRUE) or
## vcov(....) if you need it##
## $TIAM
## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: cbind(survs, deaths) ~ site + pesticide + exclosure + census +
## A.con_curt + A.het_curt + dens + (1 | quad.unique)
## Data: spdat
## Control:
## glmerControl(optimizer = "optimx", optCtrl = list(method = c("bobyqa",
## "Nelder-Mead")))
##
## AIC BIC logLik deviance df.resid
## 1992.2 2055.4 -982.1 1964.2 659
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.4000 -0.6972 0.0891 0.8041 2.1376
##
## Random effects:
## Groups Name Variance Std.Dev.
## quad.unique (Intercept) 0.2655 0.5153
## Number of obs: 673, groups: quad.unique, 233
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.66161 1.02028 1.629 0.10340
## site2 -0.16619 0.14091 -1.179 0.23826
## site3 0.23276 0.14022 1.660 0.09693 .
## pesticideF 0.38856 0.14601 2.661 0.00779 **
## pesticideI 0.14049 0.15422 0.911 0.36231
## pesticideC -0.18497 0.15618 -1.184 0.23628
## exclosure1 0.16308 0.11200 1.456 0.14536
## census16sp -1.40203 0.46003 -3.048 0.00231 **
## census17fa -0.25289 0.09225 -2.741 0.00612 **
## census17sp -1.17234 0.10043 -11.673 < 2e-16 ***
## A.con_curt -0.20009 0.04056 -4.933 8.1e-07 ***
## A.het_curt 0.08584 0.05139 1.670 0.09483 .
## dens -0.26475 0.08656 -3.059 0.00222 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1##
## Correlation matrix not shown by default, as p = 13 > 12.
## Use print(...., correlation=TRUE) or
## vcov(....) if you need itlapply(nddmods, function(x) anova(x))## $ACBA
## Analysis of Variance Table
## Df Sum Sq Mean Sq F value
## site 2 2.6593 1.3296 1.3296
## pesticide 3 0.6135 0.2045 0.2045
## exclosure 1 1.0220 1.0220 1.0220
## census 3 14.4767 4.8256 4.8256
## A.con_curt 1 0.0015 0.0015 0.0015
## A.het_curt 1 0.0424 0.0424 0.0424
## dens 1 2.9350 2.9350 2.9350
##
## $ACPS
## Analysis of Variance Table
## Df Sum Sq Mean Sq F value
## site 2 0.6643 0.3322 0.3322
## pesticide 3 0.3228 0.1076 0.1076
## exclosure 1 0.8021 0.8021 0.8021
## census 3 15.2684 5.0895 5.0895
## A.con_curt 1 0.4563 0.4563 0.4563
## A.het_curt 1 0.0728 0.0728 0.0728
## dens 1 1.0919 1.0919 1.0919
##
## $ACTE
## Analysis of Variance Table
## Df Sum Sq Mean Sq F value
## site 1 0.2896 0.2896 0.2896
## pesticide 3 2.5762 0.8587 0.8587
## exclosure 1 0.1924 0.1924 0.1924
## census 3 8.5964 2.8655 2.8655
## A.con_curt 1 3.5758 3.5758 3.5758
## A.het_curt 1 5.2620 5.2620 5.2620
## dens 1 0.0456 0.0456 0.0456
##
## $FRMA
## Analysis of Variance Table
## Df Sum Sq Mean Sq F value
## site 2 10.127 5.063 5.0634
## pesticide 3 4.128 1.376 1.3759
## exclosure 1 4.588 4.588 4.5884
## census 3 139.199 46.400 46.3997
## A.con_curt 1 0.326 0.326 0.3261
## A.het_curt 1 3.534 3.534 3.5341
## dens 1 2.478 2.478 2.4781
##
## $PIKO
## Analysis of Variance Table
## Df Sum Sq Mean Sq F value
## site 2 0.0797 0.0399 0.0399
## pesticide 3 9.7310 3.2437 3.2437
## exclosure 1 0.0473 0.0473 0.0473
## census 3 3.1894 1.0631 1.0631
## A.con_curt 1 2.1612 2.1612 2.1612
## A.het_curt 1 0.0727 0.0727 0.0727
## dens 1 0.0192 0.0192 0.0192
##
## $TIAM
## Analysis of Variance Table
## Df Sum Sq Mean Sq F value
## site 2 10.486 5.243 5.2429
## pesticide 3 9.284 3.095 3.0947
## exclosure 1 3.021 3.021 3.0215
## census 3 143.101 47.700 47.7003
## A.con_curt 1 33.836 33.836 33.8358
## A.het_curt 1 2.565 2.565 2.5647
## dens 1 9.901 9.901 9.9010##---------------- Survival of other tree seedlings-------------------##
# sp.sel2 <- c('ACBA', 'ACPS', 'ACTE', 'FRMA', 'PIKO', 'TIAM')
## other tree species
pestdat.tree_adult1 <- pestdat.ag_adult[pestdat.ag_adult$growth.form == 'tree',][!(pestdat.ag_adult[pestdat.ag_adult$growth.form == 'tree',]$sp. %in% sp.sel2),]
summary(pestdat.tree_adult1)## site quadrat sp. pesticide exclosure growth.form
## 1:139 Min. :101.0 ACMO :139 W: 83 0:172 shrub: 0
## 2:198 1st Qu.:206.0 ACMA : 80 F:126 1:203 tree :375
## 3: 38 Median :305.0 QUMO : 52 I: 97
## Mean :315.8 ABNE : 46 C: 69
## 3rd Qu.:410.0 MAAM : 12
## Max. :510.0 ACTR : 11
## (Other): 35
## census deaths total survs
## 16fa: 76 Min. :0.0000 Min. : 1.000 Min. :0.000
## 16sp: 37 1st Qu.:0.0000 1st Qu.: 1.000 1st Qu.:1.000
## 17fa:189 Median :0.0000 Median : 1.000 Median :1.000
## 17sp: 73 Mean :0.3387 Mean : 1.555 Mean :1.216
## 3rd Qu.:0.0000 3rd Qu.: 1.000 3rd Qu.:1.000
## Max. :6.0000 Max. :12.000 Max. :9.000
##
## quad.unique plot quad.sp A.sum
## 2_1_209: 15 1_0: 50 1_0_105_ACTR: 4 Min. :4071
## 2_0_206: 8 1_1: 89 1_0_406_ACMA: 4 1st Qu.:4906
## 1_1_407: 7 2_0:103 1_1_305_ACMA: 4 Median :5403
## 1_1_410: 7 2_1: 95 1_1_308_ACMO: 4 Mean :5562
## 2_0_204: 7 3_0: 19 1_1_406_ACMO: 4 3rd Qu.:6041
## 2_1_206: 7 3_1: 19 1_1_407_ACMO: 4 Max. :8723
## (Other):324 (Other) :351
## A.con A.het con.dens A.con_curt
## Min. : 0.000 Min. :2865 Min. : 0.000 Min. : 0.000
## 1st Qu.: 8.913 1st Qu.:4499 1st Qu.: 1.000 1st Qu.: 2.073
## Median : 206.523 Median :5105 Median : 6.000 Median : 5.911
## Mean : 344.675 Mean :5218 Mean : 7.691 Mean : 5.172
## 3rd Qu.: 358.975 3rd Qu.:5694 3rd Qu.: 9.000 3rd Qu.: 7.107
## Max. :3114.708 Max. :8560 Max. :46.000 Max. :14.604
##
## A.het_curt A.sum_curt con.dens_curt dens
## Min. :14.20 Min. :15.97 Min. :0.000 Min. :-0.3037
## 1st Qu.:16.51 1st Qu.:16.99 1st Qu.:1.000 1st Qu.:-0.3037
## Median :17.22 Median :17.55 Median :1.817 Median :-0.3037
## Mean :17.27 Mean :17.67 Mean :1.510 Mean :-0.2098
## 3rd Qu.:17.86 3rd Qu.:18.21 3rd Qu.:2.080 3rd Qu.:-0.3037
## Max. :20.46 Max. :20.59 Max. :3.583 Max. : 1.5574
## mod.s.t_adult.r1 <- glmer(cbind(survs, deaths) ~ site + pesticide + exclosure + census +
A.con_curt + A.het_curt + (1|quad.unique) + (1|sp.),
data=pestdat.tree_adult1, family=binomial,
control=glmerControl(optimizer="optimx",
optCtrl=list(method=c('bobyqa', 'Nelder-Mead'))))
mod.s.t_adult.r2 <- glmer(cbind(survs, deaths) ~ pesticide + exclosure + census +
A.con_curt + A.het_curt + (1|plot/quad.unique) + (1|sp.),
data=pestdat.tree_adult1, family=binomial,
control=glmerControl(optimizer="optimx",
optCtrl=list(method=c('bobyqa', 'Nelder-Mead'))))## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control
## $checkConv, : Model failed to converge with max|grad| = 0.046875 (tol =
## 0.001, component 1)## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, : Model is nearly unidentifiable: very large eigenvalue
## - Rescale variables?anova(mod.s.t_adult.r1, mod.s.t_adult.r2)## Data: pestdat.tree_adult1
## Models:
## mod.s.t_adult.r2: cbind(survs, deaths) ~ pesticide + exclosure + census + A.con_curt +
## mod.s.t_adult.r2: A.het_curt + (1 | plot/quad.unique) + (1 | sp.)
## mod.s.t_adult.r1: cbind(survs, deaths) ~ site + pesticide + exclosure + census +
## mod.s.t_adult.r1: A.con_curt + A.het_curt + (1 | quad.unique) + (1 | sp.)
## Df AIC BIC logLik deviance Chisq Chi Df
## mod.s.t_adult.r2 13 421.61 472.66 -197.81 395.61
## mod.s.t_adult.r1 14 422.78 477.76 -197.39 394.78 0.8275 1
## Pr(>Chisq)
## mod.s.t_adult.r2
## mod.s.t_adult.r1 0.363# first, check the random effects
sjp.lmer(mod.s.t_adult.r2, type='re', sort.est="sort.all", p.kr = FALSE)## Sorting each group of random intercept ('sort.all') is not possible when 'facet.grid = TRUE'.
## Plotting random effects...## Sorting each group of random intercept ('sort.all') is not possible when 'facet.grid = TRUE'.## Plotting random effects...
## Sorting each group of random intercept ('sort.all') is not possible when 'facet.grid = TRUE'.
## Plotting random effects...
## setting plot as random effects doesn't take any variance, the first model with site as fixed effect woulb be appropriate
mod.s.t_adult.r.diags <- augment(mod.s.t_adult.r1)
qplot(data=mod.s.t_adult.r.diags, x=.fitted, y=.resid, geom=c('point', 'smooth')) +
geom_hline(yintercept=0, linetype='dotted')## `geom_smooth()` using method = 'loess'
qplot(data=mod.s.t_adult.r.diags, x=pesticide, y=.wtres, geom="boxplot", colour=exclosure) 
## plot fixed and random effects
sjp.lmer(mod.s.t_adult.r1, type='fe', p.kr = FALSE)
sjp.lmer(mod.s.t_adult.r1, type='re', sort.est="sort.all", p.kr = FALSE)## Sorting each group of random intercept ('sort.all') is not possible when 'facet.grid = TRUE'.
## Plotting random effects...## Sorting each group of random intercept ('sort.all') is not possible when 'facet.grid = TRUE'.## Plotting random effects...
qqPlot(resid(mod.s.t_adult.r1))
summary(mod.s.t_adult.r1, correlation=FALSE)## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: cbind(survs, deaths) ~ site + pesticide + exclosure + census +
## A.con_curt + A.het_curt + (1 | quad.unique) + (1 | sp.)
## Data: pestdat.tree_adult1
## Control:
## glmerControl(optimizer = "optimx", optCtrl = list(method = c("bobyqa",
## "Nelder-Mead")))
##
## AIC BIC logLik deviance df.resid
## 422.8 477.8 -197.4 394.8 361
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.7012 0.2109 0.3113 0.4231 1.4231
##
## Random effects:
## Groups Name Variance Std.Dev.
## quad.unique (Intercept) 0.4184 0.6469
## sp. (Intercept) 0.9803 0.9901
## Number of obs: 375, groups: quad.unique, 138; sp., 14
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 4.38646 2.71314 1.617 0.10593
## site1 -0.05411 0.27860 -0.194 0.84601
## site2 -0.23912 0.28191 -0.848 0.39632
## pesticideF 0.27889 0.40408 0.690 0.49007
## pesticideI 0.09446 0.39454 0.239 0.81079
## pesticideC -0.24100 0.41618 -0.579 0.56253
## exclosure1 -0.31163 0.30125 -1.034 0.30091
## census16sp -1.64748 0.56673 -2.907 0.00365 **
## census17fa -0.28476 0.50103 -0.568 0.56981
## census17sp -0.97794 0.51302 -1.906 0.05662 .
## A.con_curt -0.07461 0.07035 -1.060 0.28890
## A.het_curt -0.08986 0.14921 -0.602 0.54702
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# no significant effect for the rest tree seedlingsOverall, fungicide slightly higher survival of tree seedlings pooled. Effect of conspecific neighboring trees is negative while effect of heterospecific neighboring trees is positive. Seedling density shows negative density dependence. When tested at species level. Results are consistent with that at community-level.Fungicide heigher survival of two out of nine abundant species (PIKO and TIAM). Conspecific neighboring trees have negative effect on seedling survival of TIAM. Heterospecific neighboring trees have positive effect on ACTE and FRAM. Negative density shows in FRMA and TIAM. Fungicide heigher suvival when the rest of tree species pooled (i.e., not including the common species).
Analysis of Tree Recruitment
Tree recruitments mainly occur in the first census (spring), we only use two spring censuses (16sp and 17sp) in testing recruitment.
###-------------------- recruitment annlysis -----------------------------###
pest.rect.ag_adult <- read.csv("data/pest.rect.ag_adult.csv")
pest.rect.ag_adult$pesticide <- factor(pest.rect.ag_adult$pesticide, levels=c('W', 'F', 'I', 'C'))
pest.rect.ag_adult$census <- factor(pest.rect.ag_adult$census)
pest.rect.ag_adult$exclosure <- factor(pest.rect.ag_adult$exclosure)
pest.rect.ag_adult$site <- factor(pest.rect.ag_adult$site)
summary(pest.rect.ag_adult)## X site quadrat sp. pesticide exclosure
## Min. : 1 1:497 Min. :101.0 TIAM :410 W:352 0:715
## 1st Qu.: 357 2:564 1st Qu.:204.0 FRMA :283 F:401 1:710
## Median : 713 3:364 Median :304.0 ACPS :198 I:357
## Mean : 713 Mean :304.9 ACBA :126 C:315
## 3rd Qu.:1069 3rd Qu.:408.0 ACMO : 70
## Max. :1425 Max. :510.0 ACTE : 64
## (Other):274
## growth.form census recruits quad.unique plot
## shrub: 113 16sp:747 Min. : 1.000 2_1_302: 13 1_0:228
## tree :1312 17sp:678 1st Qu.: 1.000 2_0_504: 12 1_1:269
## Median : 2.000 2_0_505: 12 2_0:303
## Mean : 3.331 1_1_410: 11 2_1:261
## 3rd Qu.: 5.000 2_0_310: 11 3_0:184
## Max. :21.000 1_0_410: 10 3_1:180
## (Other):1356
## quad.sp A.sum A.con A.het
## 1_0_101_ACPS: 2 Min. : 3990 Min. : 0.00 Min. : 2374
## 1_0_101_TIAM: 2 1st Qu.: 4991 1st Qu.: 75.49 1st Qu.: 4039
## 1_0_102_FRMA: 2 Median : 5463 Median : 398.78 Median : 4877
## 1_0_102_TIAM: 2 Mean : 5639 Mean : 742.11 Mean : 4897
## 1_0_103_ACMO: 2 3rd Qu.: 6084 3rd Qu.:1248.58 3rd Qu.: 5573
## 1_0_103_TIAM: 2 Max. :10411 Max. :5308.94 Max. :10411
## (Other) :1413
## con.dens A.con_curt A.het_curt A.sum_curt
## Min. : 0.00 Min. : 0.000 Min. :13.34 Min. :15.86
## 1st Qu.: 1.00 1st Qu.: 4.226 1st Qu.:15.93 1st Qu.:17.09
## Median : 8.00 Median : 7.361 Median :16.96 Median :17.61
## Mean :10.22 Mean : 7.087 Mean :16.87 Mean :17.74
## 3rd Qu.:14.00 3rd Qu.:10.768 3rd Qu.:17.73 3rd Qu.:18.26
## Max. :46.00 Max. :17.445 Max. :21.84 Max. :21.84
##
## con.dens_curt
## Min. :0.000
## 1st Qu.:1.000
## Median :2.000
## Mean :1.766
## 3rd Qu.:2.410
## Max. :3.583
## ###------------- Recruitment models ------------------###
pest.rect.ag_adult$site <- factor(pest.rect.ag_adult$site)
contrasts(pest.rect.ag_adult$site) <- "contr.sum"
## Only tree species
pest.t.rect_adult <- pest.rect.ag_adult[pest.rect.ag_adult$growth.form == 'tree',]
mod.pest.recr1 <- glmer(recruits ~ site + pesticide + exclosure + census + A.con_curt + A.het_curt +
(1|quad.unique) + (1|sp.),
data=pest.t.rect_adult, family=poisson,
control=glmerControl(optimizer="optimx",
optCtrl=list(method=c('bobyqa', 'Nelder-Mead'))))
sjp.lmer(mod.pest.recr1, type='fe')
sjp.lmer(mod.pest.recr1, type='re', sort.est="sort.all" )## Sorting each group of random intercept ('sort.all') is not possible when 'facet.grid = TRUE'.## Plotting random effects...## Sorting each group of random intercept ('sort.all') is not possible when 'facet.grid = TRUE'.## Plotting random effects...
## Random effects of FRMA and TIAM go extramely, try to put group them against other species as the fixed effect
pest.t.rect_adult$sp.group <- ifelse(pest.t.rect_adult$sp.=='FRMA' | pest.t.rect_adult$sp.=='TIAM', 1, 0)
pest.t.rect_adult$sp.group <- factor(pest.t.rect_adult$sp.group)
contrasts(pest.t.rect_adult$site) <- "contr.sum"
summary(pest.t.rect_adult)## X site quadrat sp. pesticide exclosure
## Min. : 1.0 1:452 Min. :101 TIAM :410 W:325 0:657
## 1st Qu.: 330.8 2:516 1st Qu.:204 FRMA :283 F:368 1:655
## Median : 658.5 3:344 Median :304 ACPS :198 I:327
## Mean : 664.2 Mean :305 ACBA :126 C:292
## 3rd Qu.: 986.2 3rd Qu.:408 ACMO : 70
## Max. :1425.0 Max. :510 ACTE : 64
## (Other):161
## growth.form census recruits quad.unique plot
## shrub: 0 16sp:689 Min. : 1.000 2_1_302: 12 1_0:208
## tree :1312 17sp:623 1st Qu.: 1.000 2_0_504: 11 1_1:244
## Median : 2.000 1_1_207: 10 2_0:276
## Mean : 3.528 1_1_503: 10 2_1:240
## 3rd Qu.: 5.000 2_0_505: 10 3_0:173
## Max. :21.000 2_1_209: 10 3_1:171
## (Other):1249
## quad.sp A.sum A.con A.het
## 1_0_101_ACPS: 2 Min. : 3990 Min. : 0.0 Min. : 2374
## 1_0_101_TIAM: 2 1st Qu.: 4991 1st Qu.: 159.0 1st Qu.: 3942
## 1_0_102_FRMA: 2 Median : 5463 Median : 507.2 Median : 4805
## 1_0_102_TIAM: 2 Mean : 5636 Mean : 806.0 Mean : 4830
## 1_0_103_ACMO: 2 3rd Qu.: 6086 3rd Qu.:1345.7 3rd Qu.: 5496
## 1_0_103_TIAM: 2 Max. :10411 Max. :5308.9 Max. :10186
## (Other) :1300
## con.dens A.con_curt A.het_curt A.sum_curt
## Min. : 0.0 Min. : 0.000 Min. :13.34 Min. :15.86
## 1st Qu.: 3.0 1st Qu.: 5.417 1st Qu.:15.80 1st Qu.:17.09
## Median : 9.0 Median : 7.975 Median :16.87 Median :17.61
## Mean :11.1 Mean : 7.692 Mean :16.79 Mean :17.74
## 3rd Qu.:14.0 3rd Qu.:11.040 3rd Qu.:17.65 3rd Qu.:18.26
## Max. :46.0 Max. :17.445 Max. :21.68 Max. :21.84
##
## con.dens_curt sp.group
## Min. :0.000 0:619
## 1st Qu.:1.442 1:693
## Median :2.080
## Mean :1.915
## 3rd Qu.:2.410
## Max. :3.583
## mod.pest.recr2 <- glmer(recruits ~ site + pesticide + exclosure + census +
A.con_curt + A.het_curt + sp.group +
(1|quad.unique) + (1|sp.),
data=pest.t.rect_adult, family=poisson,
control=glmerControl(optimizer="optimx",
optCtrl=list(method=c('bobyqa', 'Nelder-Mead'))))## Warning in optimx.check(par, optcfg$ufn, optcfg$ugr, optcfg$uhess, lower, : Parameters or bounds appear to have different scalings.
## This can cause poor performance in optimization.
## It is important for derivative free methods like BOBYQA, UOBYQA, NEWUOA.mod.pest.recr3 <- glmer(recruits ~ pesticide + exclosure + census +
A.con_curt + A.het_curt + sp.group +
(1|plot/quadrat) + (1|sp.),
data=pest.t.rect_adult, family=poisson,
control=glmerControl(optimizer="optimx",
optCtrl=list(method=c('bobyqa', 'Nelder-Mead'))))
anova(mod.pest.recr1, mod.pest.recr2, mod.pest.recr3)## Data: pest.t.rect_adult
## Models:
## mod.pest.recr1: recruits ~ site + pesticide + exclosure + census + A.con_curt +
## mod.pest.recr1: A.het_curt + (1 | quad.unique) + (1 | sp.)
## mod.pest.recr3: recruits ~ pesticide + exclosure + census + A.con_curt + A.het_curt +
## mod.pest.recr3: sp.group + (1 | plot/quadrat) + (1 | sp.)
## mod.pest.recr2: recruits ~ site + pesticide + exclosure + census + A.con_curt +
## mod.pest.recr2: A.het_curt + sp.group + (1 | quad.unique) + (1 | sp.)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## mod.pest.recr1 12 5921.6 5983.8 -2948.8 5897.6
## mod.pest.recr3 12 5918.0 5980.1 -2947.0 5894.0 3.6973 0 < 2e-16
## mod.pest.recr2 13 5914.3 5981.6 -2944.1 5888.3 5.6902 1 0.01706
##
## mod.pest.recr1
## mod.pest.recr3 ***
## mod.pest.recr2 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1mod.pest.recr.diags2 <- augment(mod.pest.recr2)
qplot(data=mod.pest.recr.diags2, x=.fitted, y=.resid, geom=c('point', 'smooth')) +
geom_hline(yintercept=0, linetype='dotted')## `geom_smooth()` using method = 'gam'
qplot(data = mod.pest.recr.diags2, x= pesticide, y=.wtres, geom="boxplot", colour=exclosure) 
qqPlot(resid(mod.pest.recr2))
sjp.lmer(mod.pest.recr2, type='fe')
sjp.lmer(mod.pest.recr2, type='re', sort.est="sort.all" )## Sorting each group of random intercept ('sort.all') is not possible when 'facet.grid = TRUE'.
## Plotting random effects...## Sorting each group of random intercept ('sort.all') is not possible when 'facet.grid = TRUE'.## Plotting random effects...
summary(mod.pest.recr2, correlation=FALSE)## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: poisson ( log )
## Formula: recruits ~ site + pesticide + exclosure + census + A.con_curt +
## A.het_curt + sp.group + (1 | quad.unique) + (1 | sp.)
## Data: pest.t.rect_adult
## Control:
## glmerControl(optimizer = "optimx", optCtrl = list(method = c("bobyqa",
## "Nelder-Mead")))
##
## AIC BIC logLik deviance df.resid
## 5914.3 5981.6 -2944.1 5888.3 1299
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.6577 -0.8584 -0.3143 0.4867 7.4046
##
## Random effects:
## Groups Name Variance Std.Dev.
## quad.unique (Intercept) 0.06006 0.2451
## sp. (Intercept) 0.08898 0.2983
## Number of obs: 1312, groups: quad.unique, 238; sp., 17
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.647454 0.385920 1.678 0.093407 .
## site1 -0.038577 0.032841 -1.175 0.240136
## site2 0.212409 0.031561 6.730 1.7e-11 ***
## pesticideF 0.204853 0.060803 3.369 0.000754 ***
## pesticideI 0.001089 0.063150 0.017 0.986246
## pesticideC -0.088692 0.064183 -1.382 0.167014
## exclosure1 0.002293 0.045819 0.050 0.960095
## census17sp -0.397786 0.033850 -11.751 < 2e-16 ***
## A.con_curt 0.027352 0.007757 3.526 0.000422 ***
## A.het_curt -0.011216 0.020534 -0.546 0.584906
## sp.group1 0.889435 0.238743 3.725 0.000195 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## convergence code: 0
## Parameters or bounds appear to have different scalings.
## This can cause poor performance in optimization.
## It is important for derivative free methods like BOBYQA, UOBYQA, NEWUOA.##The latter model considering fits better, but results do not's change, so it's ok to use the previous model
###---------------- individual species analysis ----------------------------###
## reps (lets say 10) and variation in abundance (at least >5
sp.rect.counts <- summarise(group_by(pest.t.rect_adult, sp.),
n.quadrat=length(recruits), max.rect=max(recruits), min.rect=min(recruits))
sp.rect.counts$range <- sp.rect.counts$max.rect - sp.rect.counts$min.rect
sp.rect.sel <- sp.rect.counts$sp.[sp.rect.counts$n.quadrat > 10 & sp.rect.counts$range > 5]
sp.rect.counts[sp.rect.counts$sp. %in% sp.rect.sel,]## # A tibble: 8 × 5
## sp. n.quadrat max.rect min.rect range
## <fctr> <int> <int> <int> <int>
## 1 ABNE 45 12 1 11
## 2 ACBA 126 13 1 12
## 3 ACPS 198 12 1 11
## 4 ACTE 64 7 1 6
## 5 FRMA 283 21 1 20
## 6 PIKO 43 8 1 7
## 7 QUMO 40 8 1 7
## 8 TIAM 410 21 1 20# sp.rect.sel2 <- c('ABNE', 'ACBA', 'ACPS', 'ACTE', 'FRMA', 'PIKO','QUMO', 'TIAM')
sp.rect.sel2 <- c('ACBA', 'ACPS', 'ACTE', 'FRMA', 'PIKO','QUMO', 'TIAM')
rectmods <- sapply(sp.rect.sel2, function(sp){
print(sp)
spdat <- filter(pest.t.rect_adult, sp. == sp)
spdat <- droplevels(spdat)
mod <- glmer(recruits ~ site + pesticide + exclosure + census +
A.con_curt + A.het_curt + (1|quad.unique),
data=spdat, family=poisson,
control=glmerControl(optimizer="optimx",
optCtrl=list(method=c('bobyqa', 'Nelder-Mead')))
)
}, simplify=FALSE)## [1] "ACBA"
## [1] "ACPS"## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control
## $checkConv, : Model failed to converge with max|grad| = 0.0323344 (tol =
## 0.001, component 1)## [1] "ACTE"
## [1] "FRMA"## Warning in optimx.check(par, optcfg$ufn, optcfg$ugr, optcfg$uhess, lower, : Parameters or bounds appear to have different scalings.
## This can cause poor performance in optimization.
## It is important for derivative free methods like BOBYQA, UOBYQA, NEWUOA.## [1] "PIKO"
## [1] "QUMO"
## [1] "TIAM"## diagnostic plots
indmods.rect.resids <- lapply(rectmods, augment)
indmods.rect.resids <- do.call('rbind', mapply(function(spnm, dat) {
dat$species <- spnm
return(dat)
}, spnm=as.list(names(indmods.rect.resids)), dat=indmods.rect.resids,
SIMPLIFY=FALSE))
ggplot(data=indmods.rect.resids, aes(x=.fitted, y=.wtres, colour=pesticide)) +
geom_jitter() + geom_smooth() + facet_wrap(~species, scales='free') +
geom_hline(yintercept=0, linetype='dashed')## `geom_smooth()` using method = 'loess'
lapply(rectmods, function(mod){
dotplot(ranef(mod, condVar=TRUE))})## $ACBA
## $ACBA$quad.unique
##
##
## $ACPS
## $ACPS$quad.unique
##
##
## $ACTE
## $ACTE$quad.unique
##
##
## $FRMA
## $FRMA$quad.unique
##
##
## $PIKO
## $PIKO$quad.unique
##
##
## $QUMO
## $QUMO$quad.unique
##
##
## $TIAM
## $TIAM$quad.unique
lapply(rectmods, function(x) summary(x, correlation=FALSE))## $ACBA
## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: poisson ( log )
## Formula: recruits ~ site + pesticide + exclosure + census + A.con_curt +
## A.het_curt + (1 | quad.unique)
## Data: spdat
## Control:
## glmerControl(optimizer = "optimx", optCtrl = list(method = c("bobyqa",
## "Nelder-Mead")))
##
## AIC BIC logLik deviance df.resid
## 406.5 437.7 -192.2 384.5 115
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -0.8758 -0.4577 -0.2019 0.2797 2.7061
##
## Random effects:
## Groups Name Variance Std.Dev.
## quad.unique (Intercept) 0.1026 0.3203
## Number of obs: 126, groups: quad.unique, 101
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.65540 1.61951 1.022 0.3067
## site2 -0.30498 0.19797 -1.540 0.1234
## site3 -0.30455 0.24054 -1.266 0.2055
## pesticideF -0.13314 0.21391 -0.622 0.5337
## pesticideI -0.08852 0.20630 -0.429 0.6679
## pesticideC -0.42020 0.22896 -1.835 0.0665 .
## exclosure1 -0.33281 0.17976 -1.851 0.0641 .
## census17sp 0.02200 0.14881 0.148 0.8825
## A.con_curt 0.06581 0.13318 0.494 0.6212
## A.het_curt -0.04269 0.08931 -0.478 0.6326
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## $ACPS
## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: poisson ( log )
## Formula: recruits ~ site + pesticide + exclosure + census + A.con_curt +
## A.het_curt + (1 | quad.unique)
## Data: spdat
## Control:
## glmerControl(optimizer = "optimx", optCtrl = list(method = c("bobyqa",
## "Nelder-Mead")))
##
## AIC BIC logLik deviance df.resid
## 710.8 746.9 -344.4 688.8 187
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.3584 -0.5834 -0.1737 0.3254 2.7500
##
## Random effects:
## Groups Name Variance Std.Dev.
## quad.unique (Intercept) 0.08081 0.2843
## Number of obs: 198, groups: quad.unique, 157
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.547614 0.002990 183.16 <2e-16 ***
## site2 0.001170 0.002985 0.39 0.695
## site3 -0.200816 0.002985 -67.27 <2e-16 ***
## pesticideF 0.367167 0.002986 122.97 <2e-16 ***
## pesticideI 0.200777 0.002985 67.25 <2e-16 ***
## pesticideC -0.111387 0.002984 -37.32 <2e-16 ***
## exclosure1 0.406537 0.002988 136.05 <2e-16 ***
## census17sp 0.631224 0.002989 211.20 <2e-16 ***
## A.con_curt 0.028199 0.002926 9.64 <2e-16 ***
## A.het_curt -0.043298 0.002237 -19.35 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## convergence code: 0
## Model failed to converge with max|grad| = 0.0323344 (tol = 0.001, component 1)
##
##
## $ACTE
## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: poisson ( log )
## Formula: recruits ~ site + pesticide + exclosure + census + A.con_curt +
## A.het_curt + (1 | quad.unique)
## Data: spdat
## Control:
## glmerControl(optimizer = "optimx", optCtrl = list(method = c("bobyqa",
## "Nelder-Mead")))
##
## AIC BIC logLik deviance df.resid
## 195.3 216.9 -87.7 175.3 54
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -0.9890 -0.3944 -0.1977 0.1924 3.7139
##
## Random effects:
## Groups Name Variance Std.Dev.
## quad.unique (Intercept) 0 0
## Number of obs: 64, groups: quad.unique, 53
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.86641 2.09376 0.891 0.373
## site2 0.13401 0.22485 0.596 0.551
## pesticideF -0.39642 0.30706 -1.291 0.197
## pesticideI -0.17990 0.27248 -0.660 0.509
## pesticideC 0.00579 0.31724 0.018 0.985
## exclosure1 -0.19940 0.33783 -0.590 0.555
## census17sp 0.16891 0.21189 0.797 0.425
## A.con_curt 0.24979 0.16561 1.508 0.131
## A.het_curt -0.13422 0.11351 -1.183 0.237
##
## $FRMA
## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: poisson ( log )
## Formula: recruits ~ site + pesticide + exclosure + census + A.con_curt +
## A.het_curt + (1 | quad.unique)
## Data: spdat
## Control:
## glmerControl(optimizer = "optimx", optCtrl = list(method = c("bobyqa",
## "Nelder-Mead")))
##
## AIC BIC logLik deviance df.resid
## 1274.0 1314.1 -626.0 1252.0 272
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.42007 -0.51942 -0.03114 0.42931 3.05340
##
## Random effects:
## Groups Name Variance Std.Dev.
## quad.unique (Intercept) 0.1421 0.377
## Number of obs: 283, groups: quad.unique, 225
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.1996672 0.6327376 -0.316 0.7523
## site2 0.8794064 0.1028044 8.554 < 2e-16 ***
## site3 0.4617114 0.1104352 4.181 2.9e-05 ***
## pesticideF -0.0304349 0.1090314 -0.279 0.7801
## pesticideI 0.0100566 0.1095858 0.092 0.9269
## pesticideC -0.0768766 0.1105538 -0.695 0.4868
## exclosure1 0.1397609 0.0801840 1.743 0.0813 .
## census17sp -1.5931387 0.1184738 -13.447 < 2e-16 ***
## A.con_curt 0.0003858 0.0115886 0.033 0.9734
## A.het_curt 0.0738762 0.0357336 2.067 0.0387 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## convergence code: 0
## Parameters or bounds appear to have different scalings.
## This can cause poor performance in optimization.
## It is important for derivative free methods like BOBYQA, UOBYQA, NEWUOA.
##
##
## $PIKO
## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: poisson ( log )
## Formula: recruits ~ site + pesticide + exclosure + census + A.con_curt +
## A.het_curt + (1 | quad.unique)
## Data: spdat
## Control:
## glmerControl(optimizer = "optimx", optCtrl = list(method = c("bobyqa",
## "Nelder-Mead")))
##
## AIC BIC logLik deviance df.resid
## 168.8 188.1 -73.4 146.8 32
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.0326 -0.4647 -0.2246 0.1686 2.0742
##
## Random effects:
## Groups Name Variance Std.Dev.
## quad.unique (Intercept) 0.1135 0.3369
## Number of obs: 43, groups: quad.unique, 42
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 4.55983 2.90398 1.570 0.1164
## site2 -0.31798 0.33577 -0.947 0.3436
## site3 0.24589 0.35496 0.693 0.4885
## pesticideF -0.02219 0.34243 -0.065 0.9483
## pesticideI -0.40482 0.37751 -1.072 0.2836
## pesticideC -0.24694 0.41616 -0.593 0.5529
## exclosure1 -0.36151 0.26560 -1.361 0.1735
## census17sp -0.58356 0.28115 -2.076 0.0379 *
## A.con_curt -0.25522 0.18845 -1.354 0.1756
## A.het_curt -0.04435 0.09835 -0.451 0.6520
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## $QUMO
## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: poisson ( log )
## Formula: recruits ~ site + pesticide + exclosure + census + A.con_curt +
## A.het_curt + (1 | quad.unique)
## Data: spdat
## Control:
## glmerControl(optimizer = "optimx", optCtrl = list(method = c("bobyqa",
## "Nelder-Mead")))
##
## AIC BIC logLik deviance df.resid
## 143.6 162.2 -60.8 121.6 29
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.1079 -0.4019 -0.1627 0.2458 2.9716
##
## Random effects:
## Groups Name Variance Std.Dev.
## quad.unique (Intercept) 0 0
## Number of obs: 40, groups: quad.unique, 40
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -3.83219 4.02768 -0.952 0.3414
## site2 -0.40959 0.62427 -0.656 0.5118
## site3 -0.66936 0.36917 -1.813 0.0698 .
## pesticideF 0.20670 0.31899 0.648 0.5170
## pesticideI 0.44270 0.30623 1.446 0.1483
## pesticideC -0.06282 0.42210 -0.149 0.8817
## exclosure1 -0.33106 0.44759 -0.740 0.4595
## census17sp -0.16427 0.50083 -0.328 0.7429
## A.con_curt 0.14229 0.11441 1.244 0.2136
## A.het_curt 0.21224 0.21168 1.003 0.3160
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## $TIAM
## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: poisson ( log )
## Formula: recruits ~ site + pesticide + exclosure + census + A.con_curt +
## A.het_curt + (1 | quad.unique)
## Data: spdat
## Control:
## glmerControl(optimizer = "optimx", optCtrl = list(method = c("bobyqa",
## "Nelder-Mead")))
##
## AIC BIC logLik deviance df.resid
## 2086.3 2130.5 -1032.1 2064.3 399
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.4275 -0.8542 -0.2085 0.5745 5.1464
##
## Random effects:
## Groups Name Variance Std.Dev.
## quad.unique (Intercept) 0.1257 0.3545
## Number of obs: 410, groups: quad.unique, 232
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.87154 0.63050 1.382 0.16688
## site2 0.07950 0.08648 0.919 0.35796
## site3 -0.25182 0.08783 -2.867 0.00414 **
## pesticideF 0.38493 0.08891 4.329 1.5e-05 ***
## pesticideI -0.01172 0.09504 -0.123 0.90187
## pesticideC -0.09726 0.09680 -1.005 0.31503
## exclosure1 -0.18156 0.06896 -2.633 0.00847 **
## census17sp -0.43175 0.04661 -9.263 < 2e-16 ***
## A.con_curt 0.04445 0.02466 1.802 0.07147 .
## A.het_curt 0.02319 0.03190 0.727 0.46737
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1lapply(rectmods, function(x) anova(x))## $ACBA
## Analysis of Variance Table
## Df Sum Sq Mean Sq F value
## site 2 3.1207 1.5604 1.5604
## pesticide 3 3.0695 1.0232 1.0232
## exclosure 1 4.5944 4.5944 4.5944
## census 1 0.0575 0.0575 0.0575
## A.con_curt 1 0.3377 0.3377 0.3377
## A.het_curt 1 0.2235 0.2235 0.2235
##
## $ACPS
## Analysis of Variance Table
## Df Sum Sq Mean Sq F value
## site 2 2.8820 1.4410 1.4410
## pesticide 3 14.9578 4.9859 4.9859
## exclosure 1 16.0084 16.0084 16.0084
## census 1 23.6082 23.6082 23.6082
## A.con_curt 1 0.2093 0.2093 0.2093
## A.het_curt 1 0.5247 0.5247 0.5247
##
## $ACTE
## Analysis of Variance Table
## Df Sum Sq Mean Sq F value
## site 1 0.04332 0.04332 0.0433
## pesticide 3 2.14198 0.71399 0.7140
## exclosure 1 2.35134 2.35134 2.3513
## census 1 0.46046 0.46046 0.4605
## A.con_curt 1 2.44721 2.44721 2.4472
## A.het_curt 1 1.39873 1.39873 1.3987
##
## $FRMA
## Analysis of Variance Table
## Df Sum Sq Mean Sq F value
## site 2 80.363 40.182 40.1817
## pesticide 3 0.819 0.273 0.2729
## exclosure 1 3.003 3.003 3.0031
## census 1 177.048 177.048 177.0479
## A.con_curt 1 0.774 0.774 0.7744
## A.het_curt 1 4.268 4.268 4.2678
##
## $PIKO
## Analysis of Variance Table
## Df Sum Sq Mean Sq F value
## site 2 2.3499 1.1749 1.1749
## pesticide 3 3.0691 1.0230 1.0230
## exclosure 1 2.7894 2.7894 2.7894
## census 1 3.2443 3.2443 3.2443
## A.con_curt 1 1.6280 1.6280 1.6280
## A.het_curt 1 0.1998 0.1998 0.1998
##
## $QUMO
## Analysis of Variance Table
## Df Sum Sq Mean Sq F value
## site 2 3.7176 1.85881 1.8588
## pesticide 3 3.6884 1.22945 1.2295
## exclosure 1 0.0554 0.05535 0.0554
## census 1 0.0109 0.01088 0.0109
## A.con_curt 1 0.6959 0.69586 0.6959
## A.het_curt 1 1.0059 1.00590 1.0059
##
## $TIAM
## Analysis of Variance Table
## Df Sum Sq Mean Sq F value
## site 2 12.745 6.372 6.3724
## pesticide 3 29.856 9.952 9.9521
## exclosure 1 6.472 6.472 6.4720
## census 1 83.198 83.198 83.1976
## A.con_curt 1 2.744 2.744 2.7437
## A.het_curt 1 0.525 0.525 0.5248###----------------- Pooling other species -------------------------###
# sp.rect.sel2 <- c('ACBA', 'ACPS', 'ACTE', 'FRMA', 'PIKO','QUMO', 'TIAM')
pest.t.rect.other_adult <- pest.t.rect_adult[!(pest.t.rect_adult$sp. %in% sp.rect.sel2),]
summary(pest.t.rect.other_adult)## X site quadrat sp. pesticide exclosure
## Min. : 21.0 1:44 Min. :101.0 ACMO :70 W:39 0:65
## 1st Qu.: 406.8 2:92 1st Qu.:206.0 ABNE :45 F:46 1:83
## Median : 587.5 3:12 Median :302.5 ACMA :19 I:37
## Mean : 620.1 Mean :304.4 MAAM : 6 C:26
## 3rd Qu.: 792.5 3rd Qu.:408.0 ULJA : 3
## Max. :1424.0 Max. :510.0 ACTR : 1
## (Other): 4
## growth.form census recruits quad.unique plot
## shrub: 0 16sp:49 Min. : 1.000 2_0_204: 4 1_0:17
## tree :148 17sp:99 1st Qu.: 1.000 2_1_209: 4 1_1:27
## Median : 1.000 1_1_503: 3 2_0:45
## Mean : 1.743 2_0_206: 3 2_1:47
## 3rd Qu.: 2.000 2_0_504: 3 3_0: 3
## Max. :12.000 2_1_206: 3 3_1: 9
## (Other):128
## quad.sp A.sum A.con A.het
## 1_0_103_ACMO: 2 Min. :4106 Min. : 0.0 Min. :3560
## 1_1_503_ACMO: 2 1st Qu.:4949 1st Qu.: 120.3 1st Qu.:4640
## 2_0_504_ACMO: 2 Median :5397 Median : 268.4 Median :5182
## 2_1_105_ACMO: 2 Mean :5620 Mean : 294.2 Mean :5326
## 2_1_209_ULJA: 2 3rd Qu.:6249 3rd Qu.: 394.1 3rd Qu.:5922
## 2_1_302_ACMO: 2 Max. :8408 Max. :1264.4 Max. :8040
## (Other) :136
## con.dens A.con_curt A.het_curt A.sum_curt
## Min. : 0.00 Min. : 0.000 Min. :15.27 Min. :16.01
## 1st Qu.: 4.75 1st Qu.: 4.937 1st Qu.:16.68 1st Qu.:17.04
## Median : 9.00 Median : 6.450 Median :17.31 Median :17.54
## Mean :13.53 Mean : 5.834 Mean :17.40 Mean :17.73
## 3rd Qu.:20.00 3rd Qu.: 7.331 3rd Qu.:18.09 3rd Qu.:18.42
## Max. :46.00 Max. :10.813 Max. :20.03 Max. :20.33
##
## con.dens_curt sp.group
## Min. :0.000 0:148
## 1st Qu.:1.679 1: 0
## Median :2.080
## Mean :2.043
## 3rd Qu.:2.714
## Max. :3.583
## mod.pest.recr.other <- glmer(recruits ~ site + census + pesticide + exclosure + A.con_curt + A.het_curt +
(1|quad.unique) + (1|sp.),
data=pest.t.rect.other_adult, family=poisson,
control=glmerControl(optimizer="optimx",
optCtrl=list(method=c('bobyqa', 'Nelder-Mead'))))## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control
## $checkConv, : Model failed to converge with max|grad| = 0.0105948 (tol =
## 0.001, component 1)mod.pest.recr.other.diags <- augment(mod.pest.recr.other)
qplot(data=mod.pest.recr.other.diags, x=.fitted, y=.resid, geom=c('point', 'smooth')) +
geom_hline(yintercept=0, linetype='dotted')## `geom_smooth()` using method = 'loess'
qplot(data = mod.pest.recr.other.diags, x= pesticide, y=.wtres, geom="boxplot", colour=exclosure) 
qqPlot(resid(mod.pest.recr.other))
sjp.lmer(mod.pest.recr.other, type='fe')
sjp.lmer(mod.pest.recr.other, type='re', sort.est="sort.all" )## Sorting each group of random intercept ('sort.all') is not possible when 'facet.grid = TRUE'.
## Plotting random effects...## Sorting each group of random intercept ('sort.all') is not possible when 'facet.grid = TRUE'.## Plotting random effects...
summary(mod.pest.recr.other, correlation=FALSE)## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: poisson ( log )
## Formula: recruits ~ site + census + pesticide + exclosure + A.con_curt +
## A.het_curt + (1 | quad.unique) + (1 | sp.)
## Data: pest.t.rect.other_adult
## Control:
## glmerControl(optimizer = "optimx", optCtrl = list(method = c("bobyqa",
## "Nelder-Mead")))
##
## AIC BIC logLik deviance df.resid
## 459.6 495.5 -217.8 435.6 136
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.5560 -0.3930 -0.1669 0.1491 4.6519
##
## Random effects:
## Groups Name Variance Std.Dev.
## quad.unique (Intercept) 3.572e-05 0.005977
## sp. (Intercept) 5.487e-02 0.234235
## Number of obs: 148, groups: quad.unique, 103; sp., 10
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.832636 0.007448 246.06 < 2e-16 ***
## site1 -0.007413 0.007682 -0.96 0.335
## site2 0.206950 0.007441 27.81 < 2e-16 ***
## census17sp 0.241130 0.007686 31.37 < 2e-16 ***
## pesticideF -0.179413 0.007683 -23.35 < 2e-16 ***
## pesticideI 0.063588 0.007440 8.55 < 2e-16 ***
## pesticideC 0.397282 0.007444 53.37 < 2e-16 ***
## exclosure1 0.165368 0.007680 21.53 < 2e-16 ***
## A.con_curt 0.044311 0.007129 6.22 5.11e-10 ***
## A.het_curt -0.120519 0.005470 -22.03 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## convergence code: 0
## Model failed to converge with max|grad| = 0.0105948 (tol = 0.001, component 1)Overall, fungicide higher recruitment, conspecific neigboring trees had a positive effect on recruitment. Tested at species level, fungicide increased recruitment of ACPS, TIAM. Conspecific neighboring trees increased recruitment of ACPS, slightly enhanced TIAM. Heterospecific neighboring trees increased recruitment of FRMA but decreased ACPS. Exclosure enhanced recruitment of ACPS. Pesticide had no effect on recruitment of the rest of species.
Growth analysis
Growth is calculated as the absolute values from the census to the next. SND are used to calculate the growth of seedlings, and outliers are estimated and removed in the final dataset.
## read in data
pestdat.g_adult <- read.csv("data/pestdat.g_adult.csv")
summary(pestdat.g_adult)## X sp. site quadrat
## Min. : 170 FRMA :3144 Min. :1.000 Min. :101.0
## 1st Qu.: 8806 TIAM :2006 1st Qu.:1.000 1st Qu.:202.0
## Median :13889 ACPS : 571 Median :2.000 Median :304.0
## Mean :13736 EUMA : 505 Mean :1.908 Mean :305.7
## 3rd Qu.:19216 PHSC : 484 3rd Qu.:2.000 3rd Qu.:409.0
## Max. :23888 ACBA : 355 Max. :3.000 Max. :510.0
## (Other):2206
## tag pesticide exclosure growth.form census
## Min. : 1001 C:1986 Min. :0.0000 shrub:2452 15fa: 778
## 1st Qu.: 3018 F:2816 1st Qu.:0.0000 tree :6819 16fa:2214
## Median : 6005 I:2177 Median :1.0000 16sp:3088
## Mean : 5701 W:2292 Mean :0.5113 17sp:3191
## 3rd Qu.: 8023 3rd Qu.:1.0000
## Max. :10052 Max. :1.0000
##
## grow.snd quad.unique plot quad.sp
## Min. :-2.211138 2_0_510: 106 1_0:1387 2_1_310_FRMA: 52
## 1st Qu.:-0.406868 2_1_310: 96 1_1:1703 2_1_103_FRMA: 48
## Median : 0.015854 3_0_108: 87 2_0:1967 3_0_108_FRMA: 48
## Mean : 0.002379 2_1_505: 84 2_1:1981 3_1_510_FRMA: 46
## 3rd Qu.: 0.347986 2_1_103: 79 3_0:1177 2_1_506_FRMA: 44
## Max. : 2.249671 2_1_301: 79 3_1:1056 2_1_304_FRMA: 43
## (Other):8740 (Other) :8990
## A.sum A.con A.het con.dens
## Min. : 3990 Min. : 0.0 Min. : 2374 Min. : 0.000
## 1st Qu.: 5026 1st Qu.: 0.0 1st Qu.: 4296 1st Qu.: 0.000
## Median : 5554 Median : 276.8 Median : 5097 Median : 3.000
## Mean : 5719 Mean : 593.3 Mean : 5125 Mean : 6.032
## 3rd Qu.: 6213 3rd Qu.: 977.7 3rd Qu.: 5791 3rd Qu.:10.000
## Max. :10411 Max. :5308.9 Max. :10411 Max. :46.000
##
## A.con_curt A.het_curt A.sum_curt con.dens_curt
## Min. : 0.000 Min. :13.34 Min. :15.86 Min. :0.000
## 1st Qu.: 0.000 1st Qu.:16.26 1st Qu.:17.13 1st Qu.:0.000
## Median : 6.517 Median :17.21 Median :17.71 Median :1.442
## Mean : 5.734 Mean :17.13 Mean :17.82 Mean :1.239
## 3rd Qu.: 9.925 3rd Qu.:17.96 3rd Qu.:18.38 3rd Qu.:2.154
## Max. :17.445 Max. :21.84 Max. :21.84 Max. :3.583
##
## season year
## fall :2992 Min. :1.000
## spring:6279 1st Qu.:1.000
## Median :2.000
## Mean :1.583
## 3rd Qu.:2.000
## Max. :2.000
## ###----------------------- Growth models ------------------------------###
pestdat.g_adult$pesticide <- factor(pestdat.g_adult$pesticide, levels=c('W', 'F', 'I', 'C'))
pestdat.g_adult$census <- factor(pestdat.g_adult$census)
pestdat.g_adult$exclosure <- factor(pestdat.g_adult$exclosure)
pestdat.g_adult$site <- factor(pestdat.g_adult$site)
contrasts(pestdat.g_adult$site) <- "contr.sum"
## tree
mod.pest.tree.grow1 <- lmer(grow.snd ~ site + pesticide + exclosure + census + A.con_curt + A.het_curt +
(1|quad.unique) + (1|sp.),
data=subset(pestdat.g_adult, growth.form == 'tree'))
mod.pest.tree.grow2 <- lmer(grow.snd ~ site + pesticide + exclosure + census + A.con_curt + A.het_curt +
(1|quad.unique) + (1|sp.) + (1|tag),
data=subset(pestdat.g_adult, growth.form == 'tree'))
mod.pest.tree.grow3 <- lmer(grow.snd ~ site + pesticide + exclosure + census +
(1|quad.unique) + (1|sp.),
data=subset(pestdat.g_adult, growth.form == 'tree'))
mod.pest.tree.grow4 <- lmer(grow.snd ~ site + pesticide + exclosure +
(1|quad.unique) + (1|census) + (1|sp.),
data=subset(pestdat.g_adult, growth.form == 'tree'))
anova(mod.pest.tree.grow1, mod.pest.tree.grow2, mod.pest.tree.grow3, mod.pest.tree.grow4)## refitting model(s) with ML (instead of REML)## Data: subset(pestdat.g_adult, growth.form == "tree")
## Models:
## mod.pest.tree.grow4: grow.snd ~ site + pesticide + exclosure + (1 | quad.unique) +
## mod.pest.tree.grow4: (1 | census) + (1 | sp.)
## mod.pest.tree.grow3: grow.snd ~ site + pesticide + exclosure + census + (1 | quad.unique) +
## mod.pest.tree.grow3: (1 | sp.)
## mod.pest.tree.grow1: grow.snd ~ site + pesticide + exclosure + census + A.con_curt +
## mod.pest.tree.grow1: A.het_curt + (1 | quad.unique) + (1 | sp.)
## mod.pest.tree.grow2: grow.snd ~ site + pesticide + exclosure + census + A.con_curt +
## mod.pest.tree.grow2: A.het_curt + (1 | quad.unique) + (1 | sp.) + (1 | tag)
## Df AIC BIC logLik deviance Chisq Chi Df
## mod.pest.tree.grow4 11 13862 13937 -6919.8 13840
## mod.pest.tree.grow3 13 13843 13932 -6908.5 13817 22.676 2
## mod.pest.tree.grow1 15 13830 13933 -6900.2 13800 16.609 2
## mod.pest.tree.grow2 16 13832 13942 -6900.2 13800 0.000 1
## Pr(>Chisq)
## mod.pest.tree.grow4
## mod.pest.tree.grow3 1.191e-05 ***
## mod.pest.tree.grow1 0.0002474 ***
## mod.pest.tree.grow2 1.0000000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1mod.pest.tree.grow.diags <- augment(mod.pest.tree.grow1)
qplot(data=mod.pest.tree.grow.diags, x=.fitted, y=.resid, geom=c('point', 'smooth')) +
geom_hline(yintercept=0, linetype='dotted')## `geom_smooth()` using method = 'gam'
qplot(data = mod.pest.tree.grow.diags, x= pesticide, y=.wtres, geom="boxplot", colour=exclosure) 
qqPlot(resid(mod.pest.tree.grow1))
summary(mod.pest.tree.grow1, correlation=FALSE)## Linear mixed model fit by REML ['lmerMod']
## Formula: grow.snd ~ site + pesticide + exclosure + census + A.con_curt +
## A.het_curt + (1 | quad.unique) + (1 | sp.)
## Data: subset(pestdat.g_adult, growth.form == "tree")
##
## REML criterion at convergence: 13876.8
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.8929 -0.5715 -0.0315 0.5488 3.4367
##
## Random effects:
## Groups Name Variance Std.Dev.
## quad.unique (Intercept) 0.004908 0.07006
## sp. (Intercept) 0.004189 0.06472
## Residual 0.438786 0.66241
## Number of obs: 6819, groups: quad.unique, 238; sp., 19
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 0.165058 0.165519 0.997
## site1 0.022380 0.014858 1.506
## site2 0.020914 0.013392 1.562
## pesticideF 0.004157 0.025929 0.160
## pesticideI -0.023546 0.027367 -0.860
## pesticideC -0.045913 0.028072 -1.636
## exclosure1 0.004993 0.019784 0.252
## census16fa 0.130424 0.047378 2.753
## census16sp -0.347956 0.046609 -7.465
## census17sp -0.504236 0.046588 -10.823
## A.con_curt -0.009241 0.003210 -2.879
## A.het_curt 0.011426 0.008806 1.298###------------------- season and year effect? --------------------------###
unique(pestdat.g_adult$census)## [1] 15fa 16sp 16fa 17sp
## Levels: 15fa 16fa 16sp 17sppestdat.g_adult$season <- ifelse(pestdat.g_adult$census=="16sp" |
pestdat.g_adult$census=="17sp", "spring","fall")
pestdat.g_adult$year <- ifelse(pestdat.g_adult$census=="15fa" |
pestdat.g_adult$census=="16sp", "1","2")
str(pestdat.g_adult)## 'data.frame': 9271 obs. of 23 variables:
## $ X : int 170 174 184 248 253 273 280 293 297 298 ...
## $ sp. : Factor w/ 38 levels "ABNE","ACBA",..: 2 2 2 2 2 2 2 2 2 2 ...
## $ site : Factor w/ 3 levels "1","2","3": 2 1 1 1 1 1 1 1 1 1 ...
## ..- attr(*, "contrasts")= chr "contr.sum"
## $ quadrat : int 402 108 309 304 304 209 310 505 503 308 ...
## $ tag : int 2003 8005 9007 4002 4011 9008 10007 5002 3004 8002 ...
## $ pesticide : Factor w/ 4 levels "W","F","I","C": 1 1 1 3 3 2 2 3 2 4 ...
## $ exclosure : Factor w/ 2 levels "0","1": 2 1 2 1 1 2 2 1 1 2 ...
## $ growth.form : Factor w/ 2 levels "shrub","tree": 2 2 2 2 2 2 2 2 2 2 ...
## $ census : Factor w/ 4 levels "15fa","16fa",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ grow.snd : num 0.2773 -0.0239 -0.2248 1.9844 0.9802 ...
## $ quad.unique : Factor w/ 238 levels "1_0_101","1_0_102",..: 145 8 63 20 20 55 64 36 34 62 ...
## $ plot : Factor w/ 6 levels "1_0","1_1","2_0",..: 4 1 2 1 1 2 2 1 1 2 ...
## $ quad.sp : Factor w/ 1417 levels "1_0_101_ACKO",..: 997 54 414 135 135 361 424 230 216 405 ...
## $ A.sum : num 6306 5299 5626 4975 4975 ...
## $ A.con : num 2.89 2.79 1.85 7.86 7.86 ...
## $ A.het : num 6303 5296 5624 4967 4967 ...
## $ con.dens : int 2 1 1 1 1 1 2 2 1 1 ...
## $ A.con_curt : num 1.42 1.41 1.23 1.99 1.99 ...
## $ A.het_curt : num 18.5 17.4 17.8 17.1 17.1 ...
## $ A.sum_curt : num 18.5 17.4 17.8 17.1 17.1 ...
## $ con.dens_curt: num 1.26 1 1 1 1 ...
## $ season : chr "fall" "fall" "fall" "fall" ...
## $ year : chr "1" "1" "1" "1" ...mod.pest.tree.grow5 <- lmer(grow.snd ~ site + pesticide + exclosure + season*year + A.con_curt + A.het_curt +
(1|quad.unique) + (1|sp.),
data=subset(pestdat.g_adult, growth.form == 'tree'))
mod.pest.tree.grow6 <- lmer(grow.snd ~ site + pesticide + exclosure + season + year + A.con_curt + A.het_curt +
(1|quad.unique) + (1|sp.),
data=subset(pestdat.g_adult, growth.form == 'tree'))
anova(mod.pest.tree.grow1, mod.pest.tree.grow5, mod.pest.tree.grow6)## refitting model(s) with ML (instead of REML)## Data: subset(pestdat.g_adult, growth.form == "tree")
## Models:
## mod.pest.tree.grow6: grow.snd ~ site + pesticide + exclosure + season + year + A.con_curt +
## mod.pest.tree.grow6: A.het_curt + (1 | quad.unique) + (1 | sp.)
## mod.pest.tree.grow1: grow.snd ~ site + pesticide + exclosure + census + A.con_curt +
## mod.pest.tree.grow1: A.het_curt + (1 | quad.unique) + (1 | sp.)
## mod.pest.tree.grow5: grow.snd ~ site + pesticide + exclosure + season * year + A.con_curt +
## mod.pest.tree.grow5: A.het_curt + (1 | quad.unique) + (1 | sp.)
## Df AIC BIC logLik deviance Chisq Chi Df
## mod.pest.tree.grow6 14 13860 13955 -6915.8 13832
## mod.pest.tree.grow1 15 13830 13933 -6900.2 13800 31.225 1
## mod.pest.tree.grow5 15 13830 13933 -6900.2 13800 0.000 0
## Pr(>Chisq)
## mod.pest.tree.grow6
## mod.pest.tree.grow1 2.297e-08 ***
## mod.pest.tree.grow5 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1summary(mod.pest.tree.grow5, correlation=FALSE)## Linear mixed model fit by REML ['lmerMod']
## Formula:
## grow.snd ~ site + pesticide + exclosure + season * year + A.con_curt +
## A.het_curt + (1 | quad.unique) + (1 | sp.)
## Data: subset(pestdat.g_adult, growth.form == "tree")
##
## REML criterion at convergence: 13876.8
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.8929 -0.5715 -0.0315 0.5488 3.4367
##
## Random effects:
## Groups Name Variance Std.Dev.
## quad.unique (Intercept) 0.004908 0.07006
## sp. (Intercept) 0.004189 0.06472
## Residual 0.438786 0.66241
## Number of obs: 6819, groups: quad.unique, 238; sp., 19
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 0.165058 0.165519 0.997
## site1 0.022380 0.014858 1.506
## site2 0.020914 0.013392 1.562
## pesticideF 0.004157 0.025929 0.160
## pesticideI -0.023546 0.027367 -0.860
## pesticideC -0.045913 0.028072 -1.636
## exclosure1 0.004993 0.019784 0.252
## seasonspring -0.347956 0.046609 -7.465
## year2 0.130424 0.047378 2.753
## A.con_curt -0.009241 0.003210 -2.879
## A.het_curt 0.011426 0.008806 1.298
## seasonspring:year2 -0.286705 0.051323 -5.586## shrub
mod.pest.shrub.grow1 <- lmer(grow.snd ~ site + pesticide + exclosure + A.sum_curt + census +
(1|quad.unique) + (1|sp.),
data=subset(pestdat.g_adult, growth.form == 'shrub'))
mod.pest.shrub.grow2 <- lmer(grow.snd ~ site + pesticide + exclosure + census +
(1|quad.unique) + (1|sp.),
data=subset(pestdat.g_adult, growth.form == 'shrub'))
mod.pest.shrub.grow3 <- lmer(grow.snd ~ site + pesticide + exclosure +census +
(1|quad.unique) + (1|sp.) + (1|tag),
data=subset(pestdat.g_adult, growth.form == 'shrub'))
anova(mod.pest.shrub.grow1, mod.pest.shrub.grow2, mod.pest.shrub.grow3)## refitting model(s) with ML (instead of REML)## Data: subset(pestdat.g_adult, growth.form == "shrub")
## Models:
## mod.pest.shrub.grow2: grow.snd ~ site + pesticide + exclosure + census + (1 | quad.unique) +
## mod.pest.shrub.grow2: (1 | sp.)
## mod.pest.shrub.grow1: grow.snd ~ site + pesticide + exclosure + A.sum_curt + census +
## mod.pest.shrub.grow1: (1 | quad.unique) + (1 | sp.)
## mod.pest.shrub.grow3: grow.snd ~ site + pesticide + exclosure + census + (1 | quad.unique) +
## mod.pest.shrub.grow3: (1 | sp.) + (1 | tag)
## Df AIC BIC logLik deviance Chisq Chi Df
## mod.pest.shrub.grow2 13 5139.8 5215.3 -2556.9 5113.8
## mod.pest.shrub.grow1 14 5141.7 5222.9 -2556.8 5113.7 0.1453 1
## mod.pest.shrub.grow3 14 5141.8 5223.1 -2556.9 5113.8 0.0000 0
## Pr(>Chisq)
## mod.pest.shrub.grow2
## mod.pest.shrub.grow1 0.703
## mod.pest.shrub.grow3 1.000mod.pest.shrub.grow.diags <- augment(mod.pest.shrub.grow2)
qplot(data=mod.pest.shrub.grow.diags, x=.fitted, y=.resid, geom=c('point', 'smooth')) +
geom_hline(yintercept=0, linetype='dotted')## `geom_smooth()` using method = 'gam'
qplot(data = mod.pest.shrub.grow.diags, x= pesticide, y=.wtres, geom="boxplot", colour=exclosure) 
qqPlot(resid(mod.pest.shrub.grow2))
summary(mod.pest.shrub.grow2, correlation=FALSE)## Linear mixed model fit by REML ['lmerMod']
## Formula:
## grow.snd ~ site + pesticide + exclosure + census + (1 | quad.unique) +
## (1 | sp.)
## Data: subset(pestdat.g_adult, growth.form == "shrub")
##
## REML criterion at convergence: 5167.1
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.5074 -0.4832 -0.0120 0.5113 3.4326
##
## Random effects:
## Groups Name Variance Std.Dev.
## quad.unique (Intercept) 0.0000 0.0000
## sp. (Intercept) 0.0000 0.0000
## Residual 0.4732 0.6879
## Number of obs: 2452, groups: quad.unique, 192; sp., 19
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 0.168721 0.040591 4.157
## site1 -0.001252 0.019558 -0.064
## site2 -0.007465 0.020397 -0.366
## pesticideF -0.023955 0.038375 -0.624
## pesticideI 0.019500 0.039790 0.490
## pesticideC -0.014190 0.039680 -0.358
## exclosure1 0.014623 0.027928 0.524
## census16fa 0.057311 0.040608 1.411
## census16sp -0.272208 0.040262 -6.761
## census17sp -0.342546 0.039550 -8.661##------------------------- individual species -------------------------------##
## And now pulling the species apart.
## First need to filter down to species that can be analysed - need enough
## reps (lets say 5) and variation in abundance (at least 2x
sp.grow.counts <- summarise(group_by(pestdat.g_adult, sp.),
n.indiv=length(grow.snd), n.quadrat=length(unique(grow.snd)))
sp.grow.sel <- sp.grow.counts$sp.[sp.grow.counts$n.indiv > 50 & sp.grow.counts$n.quadrat > 10]
sp.grow.counts[sp.grow.counts$sp. %in% sp.grow.sel,]## # A tibble: 19 × 3
## sp. n.indiv n.quadrat
## <fctr> <int> <int>
## 1 ACBA 355 29
## 2 ACKO 81 23
## 3 ACMA 78 22
## 4 ACMO 142 17
## 5 ACPS 571 22
## 6 ACTE 190 20
## 7 CEOR 137 31
## 8 DEGL 150 32
## 9 EUAL 231 32
## 10 EUMA 505 36
## 11 FRMA 3144 32
## 12 LOSP 94 26
## 13 PHSC 484 45
## 14 QUMO 80 18
## 15 RIMA 80 26
## 16 RIMX 158 40
## 17 SPCH 206 42
## 18 SYRE 167 28
## 19 TIAM 2006 24sp.grow.sel2 <- c('ACBA','ACKO','ACMA','ACMO','ACPS','ACTE','CEOR','DEGL','EUAL','EUMA',
'FRMA','LOSP','PHSC','QUMO','RIMX','SPCH','SYRE', 'TIAM')
growmods <- sapply(sp.grow.sel2, function(sp){
print(sp)
spdat <- filter(pestdat.g_adult, sp. == sp)
spdat <- droplevels(spdat)
mod <- lmer(grow.snd ~ site + pesticide + exclosure +
A.con_curt + A.het_curt + census + (1|quad.unique),
data=spdat)
}, simplify=FALSE)## [1] "ACBA"
## [1] "ACKO"## fixed-effect model matrix is rank deficient so dropping 2 columns / coefficients## [1] "ACMA"
## [1] "ACMO"
## [1] "ACPS"
## [1] "ACTE"
## [1] "CEOR"## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient## [1] "DEGL"## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient## [1] "EUAL"## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient## [1] "EUMA"## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient## [1] "FRMA"
## [1] "LOSP"## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient## [1] "PHSC"## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient## [1] "QUMO"
## [1] "RIMX"## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient## [1] "SPCH"## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient## [1] "SYRE"
## [1] "TIAM"## diagnostic plots
indmods.grow.resids <- lapply(growmods, augment)
indmods.grow.resids <- do.call('rbind', mapply(function(spnm, dat) {
dat$species <- spnm
return(dat)
}, spnm=as.list(names(indmods.grow.resids)), dat=indmods.grow.resids,
SIMPLIFY=FALSE))
ggplot(data=indmods.grow.resids, aes(x=.fitted, y=.wtres, colour=pesticide)) +
geom_jitter() + geom_smooth() + facet_wrap(~species, scales='free') +
geom_hline(yintercept=0, linetype='dashed')## `geom_smooth()` using method = 'loess'
library(lattice)
lapply(growmods, function(mod){
dotplot(ranef(mod, condVar=TRUE))})## $ACBA
## $ACBA$quad.unique
##
##
## $ACKO
## $ACKO$quad.unique
##
##
## $ACMA
## $ACMA$quad.unique
##
##
## $ACMO
## $ACMO$quad.unique
##
##
## $ACPS
## $ACPS$quad.unique
##
##
## $ACTE
## $ACTE$quad.unique
##
##
## $CEOR
## $CEOR$quad.unique
##
##
## $DEGL
## $DEGL$quad.unique
##
##
## $EUAL
## $EUAL$quad.unique
##
##
## $EUMA
## $EUMA$quad.unique
##
##
## $FRMA
## $FRMA$quad.unique
##
##
## $LOSP
## $LOSP$quad.unique
##
##
## $PHSC
## $PHSC$quad.unique
##
##
## $QUMO
## $QUMO$quad.unique
##
##
## $RIMX
## $RIMX$quad.unique
##
##
## $SPCH
## $SPCH$quad.unique
##
##
## $SYRE
## $SYRE$quad.unique
##
##
## $TIAM
## $TIAM$quad.unique
lapply(growmods, function(x) summary(x, correlation=FALSE))## $ACBA
## Linear mixed model fit by REML ['lmerMod']
## Formula:
## grow.snd ~ site + pesticide + exclosure + A.con_curt + A.het_curt +
## census + (1 | quad.unique)
## Data: spdat
##
## REML criterion at convergence: 118.4
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -4.8440 -0.3841 -0.0551 0.2699 6.0125
##
## Random effects:
## Groups Name Variance Std.Dev.
## quad.unique (Intercept) 0.001069 0.0327
## Residual 0.070515 0.2655
## Number of obs: 355, groups: quad.unique, 96
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) -0.02986 0.40147 -0.074
## site2 0.01970 0.04262 0.462
## site3 -0.02626 0.05608 -0.468
## pesticideF 0.03613 0.04486 0.805
## pesticideI 0.06464 0.04259 1.518
## pesticideC 0.03307 0.04386 0.754
## exclosure1 -0.03020 0.03619 -0.835
## A.con_curt 0.02986 0.03113 0.959
## A.het_curt 0.01474 0.02201 0.670
## census16fa -0.05926 0.07281 -0.814
## census16sp -0.23297 0.07054 -3.303
## census17sp -0.24286 0.06951 -3.494
##
## $ACKO
## Linear mixed model fit by REML ['lmerMod']
## Formula:
## grow.snd ~ site + pesticide + exclosure + A.con_curt + A.het_curt +
## census + (1 | quad.unique)
## Data: spdat
##
## REML criterion at convergence: 123.5
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.2651 -0.5792 -0.0307 0.5207 2.9520
##
## Random effects:
## Groups Name Variance Std.Dev.
## quad.unique (Intercept) 0.008899 0.09433
## Residual 0.229533 0.47910
## Number of obs: 81, groups: quad.unique, 14
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) -0.17891 1.59398 -0.112
## site2 0.20542 0.32586 0.630
## site3 0.24281 0.46849 0.518
## pesticideF -0.14537 0.18264 -0.796
## pesticideI -0.22591 0.19898 -1.135
## pesticideC 0.02469 0.19601 0.126
## A.het_curt 0.01079 0.09073 0.119
## census16fa -0.15840 0.15464 -1.024
## census16sp -0.07924 0.14989 -0.529
## census17sp -0.15407 0.14989 -1.028
## fit warnings:
## fixed-effect model matrix is rank deficient so dropping 2 columns / coefficients
##
## $ACMA
## Linear mixed model fit by REML ['lmerMod']
## Formula:
## grow.snd ~ site + pesticide + exclosure + A.con_curt + A.het_curt +
## census + (1 | quad.unique)
## Data: spdat
##
## REML criterion at convergence: 153.6
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.13201 -0.52436 0.07435 0.45412 2.60520
##
## Random effects:
## Groups Name Variance Std.Dev.
## quad.unique (Intercept) 0.0000 0.0000
## Residual 0.3677 0.6064
## Number of obs: 78, groups: quad.unique, 28
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 2.083091 1.477054 1.410
## site2 0.074927 0.169845 0.441
## site3 0.080634 0.656902 0.123
## pesticideF 0.048753 0.226693 0.215
## pesticideI 0.106156 0.238536 0.445
## pesticideC 0.208450 0.300643 0.693
## exclosure1 0.135532 0.198238 0.684
## A.con_curt 0.006269 0.071242 0.088
## A.het_curt -0.144097 0.087893 -1.640
## census16fa 0.600047 0.230200 2.607
## census16sp 0.102611 0.221310 0.464
## census17sp 0.005543 0.226436 0.024
##
## $ACMO
## Linear mixed model fit by REML ['lmerMod']
## Formula:
## grow.snd ~ site + pesticide + exclosure + A.con_curt + A.het_curt +
## census + (1 | quad.unique)
## Data: spdat
##
## REML criterion at convergence: 289.1
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.2511 -0.5725 0.0954 0.6358 2.9564
##
## Random effects:
## Groups Name Variance Std.Dev.
## quad.unique (Intercept) 0.0000 0.0000
## Residual 0.3979 0.6308
## Number of obs: 142, groups: quad.unique, 64
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 0.61605 1.28091 0.481
## site2 0.20873 0.14123 1.478
## site3 -0.13429 0.18291 -0.734
## pesticideF -0.16813 0.14258 -1.179
## pesticideI -0.17208 0.16448 -1.046
## pesticideC -0.14416 0.19728 -0.731
## exclosure1 -0.09585 0.11379 -0.842
## A.con_curt 0.05493 0.06117 0.898
## A.het_curt -0.03423 0.06156 -0.556
## census16fa -0.20644 0.39926 -0.517
## census16sp -0.34968 0.39367 -0.888
## census17sp -0.41797 0.38980 -1.072
##
## $ACPS
## Linear mixed model fit by REML ['lmerMod']
## Formula:
## grow.snd ~ site + pesticide + exclosure + A.con_curt + A.het_curt +
## census + (1 | quad.unique)
## Data: spdat
##
## REML criterion at convergence: 935.1
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.2362 -0.5511 -0.0253 0.6237 4.0505
##
## Random effects:
## Groups Name Variance Std.Dev.
## quad.unique (Intercept) 1.180e-15 3.435e-08
## Residual 2.818e-01 5.308e-01
## Number of obs: 571, groups: quad.unique, 146
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 0.65122 0.47133 1.382
## site2 0.08588 0.05462 1.572
## site3 -0.20869 0.06869 -3.038
## pesticideF 0.03595 0.05972 0.602
## pesticideI -0.02598 0.06450 -0.403
## pesticideC 0.01486 0.06975 0.213
## exclosure1 -0.00874 0.05560 -0.157
## A.con_curt 0.02879 0.02960 0.972
## A.het_curt -0.02989 0.02548 -1.173
## census16fa 0.25725 0.15374 1.673
## census16sp -0.34192 0.15282 -2.237
## census17sp -0.36838 0.14550 -2.532
##
## $ACTE
## Linear mixed model fit by REML ['lmerMod']
## Formula:
## grow.snd ~ site + pesticide + exclosure + A.con_curt + A.het_curt +
## census + (1 | quad.unique)
## Data: spdat
##
## REML criterion at convergence: 205.2
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -4.0577 -0.5037 -0.0451 0.4431 4.2983
##
## Random effects:
## Groups Name Variance Std.Dev.
## quad.unique (Intercept) 0.0000 0.0000
## Residual 0.1467 0.3831
## Number of obs: 190, groups: quad.unique, 48
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) -1.009530 0.688932 -1.465
## site2 0.002721 0.069928 0.039
## pesticideF -0.032059 0.084703 -0.378
## pesticideI 0.068001 0.075892 0.896
## pesticideC -0.021987 0.084745 -0.260
## exclosure1 0.212109 0.088453 2.398
## A.con_curt 0.070913 0.040015 1.772
## A.het_curt 0.033758 0.037922 0.890
## census16fa 0.266343 0.154524 1.724
## census16sp 0.055113 0.151960 0.363
## census17sp -0.024637 0.150177 -0.164
##
## $CEOR
## Linear mixed model fit by REML ['lmerMod']
## Formula:
## grow.snd ~ site + pesticide + exclosure + A.con_curt + A.het_curt +
## census + (1 | quad.unique)
## Data: spdat
##
## REML criterion at convergence: 304.3
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.9594 -0.3373 0.0166 0.4735 2.6755
##
## Random effects:
## Groups Name Variance Std.Dev.
## quad.unique (Intercept) 2.695e-15 5.191e-08
## Residual 4.981e-01 7.058e-01
## Number of obs: 137, groups: quad.unique, 22
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 1.382241 1.500046 0.922
## site2 0.072388 0.156682 0.462
## pesticideF 0.132996 0.209171 0.636
## pesticideI 0.278383 0.174202 1.598
## pesticideC -0.012866 0.192453 -0.067
## exclosure1 0.243905 0.440942 0.553
## A.het_curt -0.074669 0.087294 -0.855
## census16fa 0.006404 0.181177 0.035
## census16sp -0.265446 0.178624 -1.486
## census17sp -0.473348 0.175663 -2.695
## fit warnings:
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
##
## $DEGL
## Linear mixed model fit by REML ['lmerMod']
## Formula:
## grow.snd ~ site + pesticide + exclosure + A.con_curt + A.het_curt +
## census + (1 | quad.unique)
## Data: spdat
##
## REML criterion at convergence: 339.9
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.3164 -0.4745 0.0839 0.5065 2.9166
##
## Random effects:
## Groups Name Variance Std.Dev.
## quad.unique (Intercept) 0.005447 0.0738
## Residual 0.509337 0.7137
## Number of obs: 150, groups: quad.unique, 34
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) -0.81802 1.09937 -0.744
## site2 0.10103 0.18156 0.556
## site3 0.02975 0.17806 0.167
## pesticideF 0.03134 0.18186 0.172
## pesticideI -0.12858 0.17152 -0.750
## pesticideC -0.25382 0.19298 -1.315
## exclosure1 -0.13297 0.17756 -0.749
## A.het_curt 0.07245 0.06605 1.097
## census16fa 0.01018 0.16870 0.060
## census16sp -0.61525 0.16735 -3.676
## census17sp -0.59917 0.16564 -3.617
## fit warnings:
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
##
## $EUAL
## Linear mixed model fit by REML ['lmerMod']
## Formula:
## grow.snd ~ site + pesticide + exclosure + A.con_curt + A.het_curt +
## census + (1 | quad.unique)
## Data: spdat
##
## REML criterion at convergence: 540.2
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.2848 -0.5162 -0.0112 0.6156 2.4034
##
## Random effects:
## Groups Name Variance Std.Dev.
## quad.unique (Intercept) 0.0000 0.00
## Residual 0.5624 0.75
## Number of obs: 231, groups: quad.unique, 46
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 2.078004 1.226454 1.694
## site2 -0.039674 0.118845 -0.334
## site3 -0.179866 0.159746 -1.126
## pesticideF 0.001971 0.149129 0.013
## pesticideI 0.157589 0.160361 0.983
## pesticideC 0.007698 0.147744 0.052
## exclosure1 0.034351 0.116000 0.296
## A.het_curt -0.102106 0.069144 -1.477
## census16fa 0.255233 0.144946 1.761
## census16sp -0.513113 0.145435 -3.528
## census17sp -0.572619 0.140878 -4.065
## fit warnings:
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
##
## $EUMA
## Linear mixed model fit by REML ['lmerMod']
## Formula:
## grow.snd ~ site + pesticide + exclosure + A.con_curt + A.het_curt +
## census + (1 | quad.unique)
## Data: spdat
##
## REML criterion at convergence: 1065.7
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.11938 -0.48812 -0.01841 0.60278 2.85305
##
## Random effects:
## Groups Name Variance Std.Dev.
## quad.unique (Intercept) 0.000 0.0000
## Residual 0.456 0.6753
## Number of obs: 505, groups: quad.unique, 81
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 0.38955 0.51048 0.763
## site2 0.01009 0.07253 0.139
## site3 0.22049 0.11242 1.961
## pesticideF 0.06040 0.08405 0.719
## pesticideI 0.11650 0.08924 1.306
## pesticideC 0.15739 0.08716 1.806
## exclosure1 -0.05267 0.06629 -0.795
## A.het_curt -0.02556 0.02988 -0.856
## census16fa 0.26496 0.08803 3.010
## census16sp -0.08835 0.08781 -1.006
## census17sp -0.08136 0.08592 -0.947
## fit warnings:
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
##
## $FRMA
## Linear mixed model fit by REML ['lmerMod']
## Formula:
## grow.snd ~ site + pesticide + exclosure + A.con_curt + A.het_curt +
## census + (1 | quad.unique)
## Data: spdat
##
## REML criterion at convergence: 6175.9
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -4.1668 -0.6048 0.0217 0.6321 3.8548
##
## Random effects:
## Groups Name Variance Std.Dev.
## quad.unique (Intercept) 1.027e-15 3.204e-08
## Residual 4.097e-01 6.401e-01
## Number of obs: 3144, groups: quad.unique, 224
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 0.376243 0.184366 2.041
## site2 -0.080184 0.033388 -2.402
## site3 -0.086539 0.035897 -2.411
## pesticideF 0.015605 0.031668 0.493
## pesticideI -0.034717 0.032590 -1.065
## pesticideC -0.060307 0.032638 -1.848
## exclosure1 -0.026709 0.023596 -1.132
## A.con_curt -0.009061 0.003407 -2.659
## A.het_curt 0.003696 0.009837 0.376
## census16fa 0.210417 0.056090 3.751
## census16sp -0.372730 0.055211 -6.751
## census17sp -0.669800 0.056171 -11.924
##
## $LOSP
## Linear mixed model fit by REML ['lmerMod']
## Formula:
## grow.snd ~ site + pesticide + exclosure + A.con_curt + A.het_curt +
## census + (1 | quad.unique)
## Data: spdat
##
## REML criterion at convergence: 205.6
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.2979 -0.5286 0.0695 0.6064 2.1696
##
## Random effects:
## Groups Name Variance Std.Dev.
## quad.unique (Intercept) 0.0000 0.0000
## Residual 0.4736 0.6882
## Number of obs: 94, groups: quad.unique, 22
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) -0.90763 1.93517 -0.469
## site2 -0.23221 0.22873 -1.015
## site3 -0.10166 0.19955 -0.509
## pesticideF 0.20728 0.21811 0.950
## pesticideI -0.10693 0.22169 -0.482
## pesticideC -0.04189 0.20987 -0.200
## exclosure1 -0.33438 0.17138 -1.951
## A.het_curt 0.07335 0.10588 0.693
## census16fa -0.07167 0.20664 -0.347
## census16sp -0.70535 0.20854 -3.382
## census17sp -0.39325 0.20287 -1.938
## fit warnings:
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
##
## $PHSC
## Linear mixed model fit by REML ['lmerMod']
## Formula:
## grow.snd ~ site + pesticide + exclosure + A.con_curt + A.het_curt +
## census + (1 | quad.unique)
## Data: spdat
##
## REML criterion at convergence: 1044.1
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.3339 -0.5485 -0.0017 0.5342 3.0945
##
## Random effects:
## Groups Name Variance Std.Dev.
## quad.unique (Intercept) 0.0000 0.0000
## Residual 0.4778 0.6912
## Number of obs: 484, groups: quad.unique, 84
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 0.5919165 0.7138356 0.829
## site2 -0.0057196 0.0709174 -0.081
## site3 -0.1024200 0.0975137 -1.050
## pesticideF -0.0002306 0.0927566 -0.002
## pesticideI 0.0176862 0.1070802 0.165
## pesticideC 0.0467419 0.0979844 0.477
## exclosure1 0.1416345 0.0641147 2.209
## A.het_curt -0.0209633 0.0413780 -0.507
## census16fa -0.2154622 0.0910099 -2.367
## census16sp -0.3959989 0.0893055 -4.434
## census17sp -0.4719985 0.0885211 -5.332
## fit warnings:
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
##
## $QUMO
## Linear mixed model fit by REML ['lmerMod']
## Formula:
## grow.snd ~ site + pesticide + exclosure + A.con_curt + A.het_curt +
## census + (1 | quad.unique)
## Data: spdat
##
## REML criterion at convergence: 147.2
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.75492 -0.57954 -0.08754 0.45398 2.25866
##
## Random effects:
## Groups Name Variance Std.Dev.
## quad.unique (Intercept) 0.1431 0.3782
## Residual 0.2485 0.4985
## Number of obs: 80, groups: quad.unique, 38
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 2.57531 2.62537 0.981
## site2 -0.20281 0.45852 -0.442
## site3 0.05050 0.28576 0.177
## pesticideF -0.04227 0.24932 -0.170
## pesticideI 0.12223 0.25590 0.478
## pesticideC -0.12616 0.32078 -0.393
## exclosure1 0.06607 0.32081 0.206
## A.con_curt -0.08607 0.08279 -1.040
## A.het_curt -0.08855 0.13988 -0.633
## census16fa -0.34245 0.58694 -0.584
## census16sp -0.66438 0.58727 -1.131
## census17sp -0.39055 0.57485 -0.679
##
## $RIMX
## Linear mixed model fit by REML ['lmerMod']
## Formula:
## grow.snd ~ site + pesticide + exclosure + A.con_curt + A.het_curt +
## census + (1 | quad.unique)
## Data: spdat
##
## REML criterion at convergence: 294.2
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.7970 -0.3281 0.0329 0.4977 2.6292
##
## Random effects:
## Groups Name Variance Std.Dev.
## quad.unique (Intercept) 0.000 0.0000
## Residual 0.334 0.5779
## Number of obs: 158, groups: quad.unique, 35
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) -1.35783 1.11062 -1.223
## site2 0.04218 0.10337 0.408
## site3 0.07470 0.14476 0.516
## pesticideF 0.03057 0.11456 0.267
## pesticideI 0.01713 0.12897 0.133
## pesticideC -0.04483 0.18020 -0.249
## exclosure1 -0.08953 0.09861 -0.908
## A.het_curt 0.08567 0.06278 1.365
## census16fa -0.02938 0.13604 -0.216
## census16sp -0.04387 0.13508 -0.325
## census17sp -0.25488 0.13201 -1.931
## fit warnings:
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
##
## $SPCH
## Linear mixed model fit by REML ['lmerMod']
## Formula:
## grow.snd ~ site + pesticide + exclosure + A.con_curt + A.het_curt +
## census + (1 | quad.unique)
## Data: spdat
##
## REML criterion at convergence: 405.1
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.1798 -0.5181 -0.0613 0.5629 2.8514
##
## Random effects:
## Groups Name Variance Std.Dev.
## quad.unique (Intercept) 1.092e-18 1.045e-09
## Residual 3.786e-01 6.153e-01
## Number of obs: 206, groups: quad.unique, 38
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 0.315818 1.165221 0.271
## site2 -0.078148 0.133826 -0.584
## site3 -0.004012 0.112478 -0.036
## pesticideF -0.119994 0.141277 -0.849
## pesticideI -0.258568 0.117254 -2.205
## pesticideC -0.085017 0.127398 -0.667
## exclosure1 -0.180097 0.100746 -1.788
## A.het_curt 0.009809 0.064988 0.151
## census16fa -0.205046 0.125477 -1.634
## census16sp -0.114337 0.123404 -0.927
## census17sp -0.460315 0.124963 -3.684
## fit warnings:
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
##
## $SYRE
## Linear mixed model fit by REML ['lmerMod']
## Formula:
## grow.snd ~ site + pesticide + exclosure + A.con_curt + A.het_curt +
## census + (1 | quad.unique)
## Data: spdat
##
## REML criterion at convergence: 370.6
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.4802 -0.4671 0.0078 0.4670 2.8430
##
## Random effects:
## Groups Name Variance Std.Dev.
## quad.unique (Intercept) 3.554e-16 1.885e-08
## Residual 4.853e-01 6.966e-01
## Number of obs: 167, groups: quad.unique, 30
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) -1.00354 1.21507 -0.826
## site2 0.20770 0.17993 1.154
## site3 0.17573 0.20076 0.875
## pesticideF -0.01741 0.18481 -0.094
## pesticideI 0.30822 0.17605 1.751
## pesticideC -0.17026 0.18533 -0.919
## exclosure1 -0.08093 0.16180 -0.500
## A.con_curt 0.04516 0.08304 0.544
## A.het_curt 0.05423 0.06468 0.838
## census16fa 0.17167 0.15853 1.083
## census16sp -0.28124 0.15766 -1.784
## census17sp -0.28498 0.15404 -1.850
##
## $TIAM
## Linear mixed model fit by REML ['lmerMod']
## Formula:
## grow.snd ~ site + pesticide + exclosure + A.con_curt + A.het_curt +
## census + (1 | quad.unique)
## Data: spdat
##
## REML criterion at convergence: 4630.8
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.01406 -0.58696 -0.06046 0.55875 2.82722
##
## Random effects:
## Groups Name Variance Std.Dev.
## quad.unique (Intercept) 0.01234 0.1111
## Residual 0.56533 0.7519
## Number of obs: 2006, groups: quad.unique, 225
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 0.642327 0.434866 1.477
## site2 0.011795 0.051174 0.230
## site3 -0.080662 0.049428 -1.632
## pesticideF 0.029144 0.050777 0.574
## pesticideI -0.009912 0.056936 -0.174
## pesticideC 0.005283 0.059945 0.088
## exclosure1 0.015908 0.040225 0.396
## A.con_curt -0.040216 0.014855 -2.707
## A.het_curt 0.024414 0.018220 1.340
## census16fa -0.292464 0.244175 -1.198
## census16sp -0.630650 0.242611 -2.599
## census17sp -0.726529 0.243035 -2.989lapply(growmods, function(x) anova(x))## $ACBA
## Analysis of Variance Table
## Df Sum Sq Mean Sq F value
## site 2 0.19044 0.09522 1.3503
## pesticide 3 0.10601 0.03534 0.5011
## exclosure 1 0.04581 0.04581 0.6497
## A.con_curt 1 0.04287 0.04287 0.6080
## A.het_curt 1 0.02771 0.02771 0.3929
## census 3 2.44742 0.81581 11.5692
##
## $ACKO
## Analysis of Variance Table
## Df Sum Sq Mean Sq F value
## site 2 0.13388 0.066938 0.2916
## pesticide 3 0.56789 0.189298 0.8247
## A.het_curt 1 0.00612 0.006123 0.0267
## census 3 0.33188 0.110625 0.4820
##
## $ACMA
## Analysis of Variance Table
## Df Sum Sq Mean Sq F value
## site 2 0.1982 0.09908 0.2695
## pesticide 3 0.0500 0.01667 0.0453
## exclosure 1 0.1599 0.15986 0.4348
## A.con_curt 1 0.0155 0.01546 0.0421
## A.het_curt 1 0.6869 0.68691 1.8683
## census 3 4.4915 1.49716 4.0720
##
## $ACMO
## Analysis of Variance Table
## Df Sum Sq Mean Sq F value
## site 2 1.62732 0.81366 2.0450
## pesticide 3 0.75575 0.25192 0.6332
## exclosure 1 0.16894 0.16894 0.4246
## A.con_curt 1 0.26634 0.26634 0.6694
## A.het_curt 1 0.06992 0.06992 0.1757
## census 3 1.22562 0.40854 1.0268
##
## $ACPS
## Analysis of Variance Table
## Df Sum Sq Mean Sq F value
## site 2 6.1451 3.0726 10.9050
## pesticide 3 0.6778 0.2259 0.8019
## exclosure 1 0.2705 0.2705 0.9600
## A.con_curt 1 0.2524 0.2524 0.8957
## A.het_curt 1 0.3625 0.3625 1.2865
## census 3 27.7756 9.2585 32.8599
##
## $ACTE
## Analysis of Variance Table
## Df Sum Sq Mean Sq F value
## site 1 0.42257 0.42257 2.8798
## pesticide 3 0.16310 0.05437 0.3705
## exclosure 1 0.53553 0.53553 3.6496
## A.con_curt 1 0.65353 0.65353 4.4538
## A.het_curt 1 0.10551 0.10551 0.7190
## census 3 2.48124 0.82708 5.6365
##
## $CEOR
## Analysis of Variance Table
## Df Sum Sq Mean Sq F value
## site 1 0.0036 0.00358 0.0072
## pesticide 3 1.7231 0.57435 1.1531
## exclosure 1 0.0012 0.00118 0.0024
## A.het_curt 1 0.3364 0.33640 0.6754
## census 3 5.6007 1.86689 3.7480
##
## $DEGL
## Analysis of Variance Table
## Df Sum Sq Mean Sq F value
## site 2 0.8091 0.4046 0.7943
## pesticide 3 0.9130 0.3043 0.5975
## exclosure 1 0.1115 0.1115 0.2189
## A.het_curt 1 0.6542 0.6542 1.2843
## census 3 13.9627 4.6542 9.1378
##
## $EUAL
## Analysis of Variance Table
## Df Sum Sq Mean Sq F value
## site 2 1.0750 0.5375 0.9556
## pesticide 3 0.7356 0.2452 0.4360
## exclosure 1 0.0588 0.0588 0.1045
## A.het_curt 1 1.3798 1.3798 2.4533
## census 3 28.5855 9.5285 16.9417
##
## $EUMA
## Analysis of Variance Table
## Df Sum Sq Mean Sq F value
## site 2 1.3905 0.6953 1.5248
## pesticide 3 1.3571 0.4524 0.9921
## exclosure 1 0.4103 0.4103 0.8999
## A.het_curt 1 0.2873 0.2873 0.6302
## census 3 10.5258 3.5086 7.6946
##
## $FRMA
## Analysis of Variance Table
## Df Sum Sq Mean Sq F value
## site 2 4.10 2.050 5.0026
## pesticide 3 2.73 0.909 2.2180
## exclosure 1 0.79 0.787 1.9208
## A.con_curt 1 4.50 4.497 10.9773
## A.het_curt 1 0.08 0.078 0.1901
## census 3 378.44 126.147 307.9073
##
## $LOSP
## Analysis of Variance Table
## Df Sum Sq Mean Sq F value
## site 2 0.2139 0.10696 0.2259
## pesticide 3 0.4119 0.13731 0.2900
## exclosure 1 1.7564 1.75638 3.7089
## A.het_curt 1 0.2958 0.29584 0.6247
## census 3 7.1127 2.37092 5.0066
##
## $PHSC
## Analysis of Variance Table
## Df Sum Sq Mean Sq F value
## site 2 0.6653 0.3327 0.6963
## pesticide 3 0.1744 0.0581 0.1217
## exclosure 1 2.2935 2.2935 4.8003
## A.het_curt 1 0.0692 0.0692 0.1447
## census 3 15.9818 5.3273 11.1500
##
## $QUMO
## Analysis of Variance Table
## Df Sum Sq Mean Sq F value
## site 2 0.01168 0.005840 0.0235
## pesticide 3 0.17480 0.058267 0.2345
## exclosure 1 0.00529 0.005290 0.0213
## A.con_curt 1 0.14338 0.143378 0.5769
## A.het_curt 1 0.09580 0.095799 0.3855
## census 3 0.54975 0.183250 0.7374
##
## $RIMX
## Analysis of Variance Table
## Df Sum Sq Mean Sq F value
## site 2 0.33322 0.16661 0.4989
## pesticide 3 0.00321 0.00107 0.0032
## exclosure 1 0.18359 0.18359 0.5497
## A.het_curt 1 0.72180 0.72180 2.1612
## census 3 1.71309 0.57103 1.7098
##
## $SPCH
## Analysis of Variance Table
## Df Sum Sq Mean Sq F value
## site 2 0.0178 0.00891 0.0235
## pesticide 3 2.0533 0.68443 1.8078
## exclosure 1 1.1318 1.13185 2.9896
## A.het_curt 1 0.0083 0.00828 0.0219
## census 3 5.7822 1.92740 5.0910
##
## $SYRE
## Analysis of Variance Table
## Df Sum Sq Mean Sq F value
## site 2 0.0373 0.01866 0.0384
## pesticide 3 5.6666 1.88888 3.8922
## exclosure 1 0.3082 0.30825 0.6352
## A.con_curt 1 0.0592 0.05924 0.1221
## A.het_curt 1 0.3271 0.32714 0.6741
## census 3 6.3363 2.11210 4.3522
##
## $TIAM
## Analysis of Variance Table
## Df Sum Sq Mean Sq F value
## site 2 1.060 0.5301 0.9377
## pesticide 3 0.502 0.1674 0.2961
## exclosure 1 0.218 0.2176 0.3849
## A.con_curt 1 7.343 7.3426 12.9882
## A.het_curt 1 1.027 1.0271 1.8168
## census 3 51.364 17.1213 30.2855##------------------------ Other species ------------------------##
pestdat.g.other_adult <- pestdat.g_adult[!(pestdat.g_adult$sp. %in% sp.grow.sel2),]
mod.pest.other.grow <- lmer(grow.snd ~ site + pesticide + exclosure + A.con_curt + A.het_curt + census +
(1|quad.unique) + (1|sp.), data=pestdat.g.other_adult)
mod.pest.other.grow.diags <- augment(mod.pest.other.grow)
qplot(data=mod.pest.other.grow.diags, x=.fitted, y=.resid, geom=c('point', 'smooth')) +
geom_hline(yintercept=0, linetype='dotted')## `geom_smooth()` using method = 'loess'
qplot(data = mod.pest.other.grow.diags, x= pesticide, y=.wtres, geom="boxplot", colour=exclosure) 
qqPlot(resid(mod.pest.other.grow))
summary(mod.pest.other.grow, correlation=FALSE)## Linear mixed model fit by REML ['lmerMod']
## Formula:
## grow.snd ~ site + pesticide + exclosure + A.con_curt + A.het_curt +
## census + (1 | quad.unique) + (1 | sp.)
## Data: pestdat.g.other_adult
##
## REML criterion at convergence: 1202.5
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.8924 -0.5571 -0.0478 0.5790 3.2409
##
## Random effects:
## Groups Name Variance Std.Dev.
## quad.unique (Intercept) 0.0007595 0.02756
## sp. (Intercept) 0.0000000 0.00000
## Residual 0.6282488 0.79262
## Number of obs: 492, groups: quad.unique, 109; sp., 20
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 0.935784 0.585889 1.597
## site1 -0.004980 0.057059 -0.087
## site2 -0.054944 0.051073 -1.076
## pesticideF -0.098299 0.107024 -0.918
## pesticideI 0.083325 0.111394 0.748
## pesticideC 0.043735 0.109637 0.399
## exclosure1 0.348849 0.079076 4.412
## A.con_curt 0.004207 0.010270 0.410
## A.het_curt -0.053095 0.032825 -1.618
## census16fa 0.141441 0.129259 1.094
## census16sp -0.302241 0.127118 -2.378
## census17sp -0.331094 0.121582 -2.723Pesticide had no effect on growth of seedlings. Results are consistent both at community- and species-level. Exclosure increased growth of the species other than the common species.
Diversity analysis
First, trees and shrubs are pooled and estimated the Shannon’ H and Inversimpson indexes. Next, trees and shrubs are estimated seperately.
library(vegan)## Loading required package: permute## This is vegan 2.4-4##----------------------- Overall --------------------------##
divers_woodpest <- summarise(
group_by(pestdat.ag_adult, site, plot, quadrat, quad.unique, pesticide, exclosure, census),
shannon=diversity(survs), simpson=diversity(survs, index='invsimpson'), species=specnumber(survs))
adult.ngbr.inf <- read.csv("data/adult.ngbr.inf.csv")
divers_woodpest$A.sum <- adult.ngbr.inf$A.sum[match(divers_woodpest$quad.unique, adult.ngbr.inf$quad.unique)]
summary(divers_woodpest)## site plot quadrat quad.unique pesticide exclosure
## 1:317 1_0:159 Min. :101.0 1_0_101: 4 W:234 0:474
## 2:318 1_1:158 1st Qu.:202.0 1_0_102: 4 F:237 1:463
## 3:302 2_0:159 Median :304.0 1_0_103: 4 I:230
## 2_1:159 Mean :304.7 1_0_104: 4 C:236
## 3_0:156 3rd Qu.:409.0 1_0_105: 4
## 3_1:146 Max. :510.0 1_0_106: 4
## (Other):913
## census shannon simpson species
## 16fa:238 Min. :0.0000 Min. :1.000 Min. : 0.000
## 16sp:223 1st Qu.:0.7356 1st Qu.:2.000 1st Qu.: 3.000
## 17fa:238 Median :1.2033 Median :2.941 Median : 4.000
## 17sp:238 Mean :1.1393 Mean : Inf Mean : 4.291
## 3rd Qu.:1.5591 3rd Qu.:4.167 3rd Qu.: 6.000
## Max. :2.2430 Max. : Inf Max. :11.000
##
## A.sum
## Min. : 3990
## 1st Qu.: 4975
## Median : 5463
## Mean : 5618
## 3rd Qu.: 6071
## Max. :10411
## divers_woodpest$A.sum_curt <- divers_woodpest$A.sum^(1/3)
ggplot(divers_woodpest, aes(x=pesticide, y=shannon, colour = exclosure)) + geom_boxplot()
divers_woodpest$census <- factor(divers_woodpest$census, levels=c('16sp', '16fa', '17sp', '17fa'))
summarise(group_by(divers_woodpest, pesticide, exclosure, census), meanH=mean(shannon),
seH=sd(shannon)/sqrt(length(shannon))) %>%
ggplot(aes(x=pesticide, y=meanH, ymin=meanH-seH, ymax=meanH+seH, colour=census,
shape=as.factor(exclosure))) + geom_pointrange() + coord_trans(y="exp") +
labs(x="pesticide", y="Effective number of species")
### shannon diversity
mod.wood.shannon1 <- lmer(shannon ~ site + pesticide*exclosure + census+ (1|quad.unique), data=divers_woodpest)
mod.wood.shannon2 <- lmer(shannon ~ site + pesticide + exclosure + census + (1|quad.unique), data=divers_woodpest)
mod.wood.shannon3 <- lmer(shannon ~ site + pesticide + exclosure + census + A.sum_curt +
(1|quad.unique), data=divers_woodpest)
mod.wood.shannon4 <- lmer(shannon ~ pesticide + exclosure + census + (1|plot/quadrat), data=divers_woodpest)
anova(mod.wood.shannon1, mod.wood.shannon2, mod.wood.shannon3, mod.wood.shannon4)## refitting model(s) with ML (instead of REML)## Data: divers_woodpest
## Models:
## mod.wood.shannon4: shannon ~ pesticide + exclosure + census + (1 | plot/quadrat)
## mod.wood.shannon2: shannon ~ site + pesticide + exclosure + census + (1 | quad.unique)
## mod.wood.shannon3: shannon ~ site + pesticide + exclosure + census + A.sum_curt +
## mod.wood.shannon3: (1 | quad.unique)
## mod.wood.shannon1: shannon ~ site + pesticide * exclosure + census + (1 | quad.unique)
## Df AIC BIC logLik deviance Chisq Chi Df
## mod.wood.shannon4 11 518.41 571.67 -248.20 496.41
## mod.wood.shannon2 12 500.14 558.25 -238.07 476.14 20.2701 1
## mod.wood.shannon3 13 500.94 563.89 -237.47 474.94 1.1952 1
## mod.wood.shannon1 15 503.22 575.86 -236.61 473.22 1.7177 2
## Pr(>Chisq)
## mod.wood.shannon4
## mod.wood.shannon2 6.724e-06 ***
## mod.wood.shannon3 0.2743
## mod.wood.shannon1 0.4237
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1plot(mod.wood.shannon2, type=c('p', 'smooth'))
library(sjPlot)
sjp.lmer(mod.wood.shannon2, type='fe')## Computing p-values via Kenward-Roger approximation. Use `p.kr = FALSE` if computation takes too long.
sjp.lmer(mod.wood.shannon2, type='re', sort.est="sort.all" )## Sorting each group of random intercept ('sort.all') is not possible when 'facet.grid = TRUE'.## Plotting random effects...
summary(mod.wood.shannon2, correlation=FALSE)## Linear mixed model fit by REML ['lmerMod']
## Formula:
## shannon ~ site + pesticide + exclosure + census + (1 | quad.unique)
## Data: divers_woodpest
##
## REML criterion at convergence: 525.5
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.7845 -0.5281 0.0199 0.5461 3.1071
##
## Random effects:
## Groups Name Variance Std.Dev.
## quad.unique (Intercept) 0.12917 0.3594
## Residual 0.05428 0.2330
## Number of obs: 937, groups: quad.unique, 239
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 0.74860 0.05639 13.28
## site1 0.24868 0.03456 7.20
## site2 0.06529 0.03456 1.89
## pesticideF 0.06764 0.06906 0.98
## pesticideI -0.00947 0.06942 -0.14
## pesticideC -0.05423 0.06907 -0.79
## exclosure1 -0.13516 0.04895 -2.76
## census16fa 0.53316 0.02184 24.41
## census17sp 0.49461 0.02184 22.64
## census17fa 0.72701 0.02184 33.28##----------------------- More diversity indices -------------------------##
### Invsimpson
range(divers_woodpest$simpson)## [1] 1 Infdivers_woodpest1 <- subset(divers_woodpest, simpson != Inf)
mod.wood.simpson1 <- lmer(simpson ~ site + pesticide*exclosure + census + (1|quad.unique), data=divers_woodpest1)
mod.wood.simpson2 <- lmer(simpson ~ site + pesticide + exclosure + census + (1|quad.unique), data=divers_woodpest1)
mod.wood.simpson3 <- lmer(simpson ~ site + pesticide + exclosure + census + A.sum_curt +
(1|quad.unique), data=divers_woodpest1)
anova(mod.wood.simpson1, mod.wood.simpson2, mod.wood.simpson3)## refitting model(s) with ML (instead of REML)## Data: divers_woodpest1
## Models:
## mod.wood.simpson2: simpson ~ site + pesticide + exclosure + census + (1 | quad.unique)
## mod.wood.simpson3: simpson ~ site + pesticide + exclosure + census + A.sum_curt +
## mod.wood.simpson3: (1 | quad.unique)
## mod.wood.simpson1: simpson ~ site + pesticide * exclosure + census + (1 | quad.unique)
## Df AIC BIC logLik deviance Chisq Chi Df
## mod.wood.simpson2 12 2610.7 2668.7 -1293.4 2586.7
## mod.wood.simpson3 13 2612.4 2675.3 -1293.2 2586.4 0.2661 1
## mod.wood.simpson1 15 2615.6 2688.1 -1292.8 2585.6 0.8460 2
## Pr(>Chisq)
## mod.wood.simpson2
## mod.wood.simpson3 0.6059
## mod.wood.simpson1 0.6551plot(mod.wood.simpson2, type=c('p', 'smooth'))
sjp.lmer(mod.wood.simpson2, type='fe')## Computing p-values via Kenward-Roger approximation. Use `p.kr = FALSE` if computation takes too long.
sjp.lmer(mod.wood.simpson2, type='re', sort.est="sort.all" )## Sorting each group of random intercept ('sort.all') is not possible when 'facet.grid = TRUE'.
## Plotting random effects...
qqPlot(resid(mod.wood.simpson2))
summary(mod.wood.simpson2, correlation=FALSE)## Linear mixed model fit by REML ['lmerMod']
## Formula:
## simpson ~ site + pesticide + exclosure + census + (1 | quad.unique)
## Data: divers_woodpest1
##
## REML criterion at convergence: 2614.5
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.9611 -0.5662 -0.0780 0.5089 4.1866
##
## Random effects:
## Groups Name Variance Std.Dev.
## quad.unique (Intercept) 0.9858 0.9929
## Residual 0.5688 0.7542
## Number of obs: 928, groups: quad.unique, 238
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 2.324130 0.160422 14.488
## site1 0.759747 0.097315 7.807
## site2 0.004584 0.097315 0.047
## pesticideF 0.201645 0.194368 1.037
## pesticideI 0.073678 0.196057 0.376
## pesticideC -0.090029 0.194237 -0.463
## exclosure1 -0.394519 0.138010 -2.859
## census16fa 0.957426 0.071367 13.415
## census17sp 1.012545 0.071462 14.169
## census17fa 1.750211 0.071409 24.510### Pielou's evenness (J)
table(divers_woodpest$species)##
## 0 1 2 3 4 5 6 7 8 9 10 11
## 9 85 118 181 139 127 116 72 56 24 8 2# need to remove those with less than 1 species, do not test this index in trees and shrubs seperately.
divers_woodpest2 <- subset(divers_woodpest, species > 1)
divers_woodpest2$evenness <- divers_woodpest2$shannon/log(divers_woodpest2$species)
mod.wood.evenness1 <- lmer(evenness ~ site + pesticide*exclosure + census +(1|quad.unique), data=divers_woodpest2)
mod.wood.evenness2 <- lmer(evenness ~ site + pesticide + exclosure + census +(1|quad.unique), data=divers_woodpest2)
mod.wood.evenness3 <- lmer(evenness ~ site + pesticide + exclosure + census +A.sum_curt +
(1|quad.unique), data=divers_woodpest2)
anova(mod.wood.evenness1, mod.wood.evenness2, mod.wood.evenness3)## refitting model(s) with ML (instead of REML)## Data: divers_woodpest2
## Models:
## mod.wood.evenness2: evenness ~ site + pesticide + exclosure + census + (1 | quad.unique)
## mod.wood.evenness3: evenness ~ site + pesticide + exclosure + census + A.sum_curt +
## mod.wood.evenness3: (1 | quad.unique)
## mod.wood.evenness1: evenness ~ site + pesticide * exclosure + census + (1 | quad.unique)
## Df AIC BIC logLik deviance Chisq Chi Df
## mod.wood.evenness2 12 -1730.1 -1673.2 877.03 -1754.1
## mod.wood.evenness3 13 -1734.4 -1672.8 880.22 -1760.4 6.3676 1
## mod.wood.evenness1 15 -1726.3 -1655.3 878.17 -1756.3 0.0000 2
## Pr(>Chisq)
## mod.wood.evenness2
## mod.wood.evenness3 0.01162 *
## mod.wood.evenness1 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1plot(mod.wood.evenness2, type=c('p', 'smooth'))
library(sjPlot)
sjp.lmer(mod.wood.evenness2, type='fe')## Computing p-values via Kenward-Roger approximation. Use `p.kr = FALSE` if computation takes too long.
sjp.lmer(mod.wood.evenness2, type='re', sort.est="sort.all" )## Sorting each group of random intercept ('sort.all') is not possible when 'facet.grid = TRUE'.
## Plotting random effects...
qqPlot(resid(mod.wood.evenness2))
summary(mod.wood.evenness2, correlation=FALSE)## Linear mixed model fit by REML ['lmerMod']
## Formula:
## evenness ~ site + pesticide + exclosure + census + (1 | quad.unique)
## Data: divers_woodpest2
##
## REML criterion at convergence: -1676.3
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -4.4809 -0.5011 0.1263 0.5480 3.3404
##
## Random effects:
## Groups Name Variance Std.Dev.
## quad.unique (Intercept) 0.005098 0.07140
## Residual 0.004794 0.06924
## Number of obs: 843, groups: quad.unique, 234
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 0.957464 0.012776 74.94
## site1 0.035170 0.007385 4.76
## site2 -0.037207 0.007385 -5.04
## pesticideF -0.013875 0.014724 -0.94
## pesticideI -0.001292 0.014897 -0.09
## pesticideC 0.001681 0.014863 0.11
## exclosure1 -0.037541 0.010535 -3.56
## census16fa -0.099753 0.007385 -13.51
## census17sp -0.067521 0.007400 -9.12
## census17fa -0.058929 0.007393 -7.97summary(mod.wood.evenness3, correlation=FALSE)## Linear mixed model fit by REML ['lmerMod']
## Formula: evenness ~ site + pesticide + exclosure + census + A.sum_curt +
## (1 | quad.unique)
## Data: divers_woodpest2
##
## REML criterion at convergence: -1674
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -4.5331 -0.5062 0.1335 0.5453 3.3869
##
## Random effects:
## Groups Name Variance Std.Dev.
## quad.unique (Intercept) 0.004952 0.07037
## Residual 0.004794 0.06924
## Number of obs: 843, groups: quad.unique, 234
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 1.200959 0.098399 12.205
## site1 0.031432 0.007451 4.218
## site2 -0.034372 0.007389 -4.652
## pesticideF -0.013587 0.014556 -0.933
## pesticideI -0.000756 0.014728 -0.051
## pesticideC 0.004779 0.014746 0.324
## exclosure1 -0.030991 0.010736 -2.887
## census16fa -0.099919 0.007384 -13.531
## census17sp -0.067687 0.007400 -9.147
## census17fa -0.059049 0.007392 -7.988
## A.sum_curt -0.013966 0.005598 -2.495###-------------------- Trees, shrubs -------------------------###
## Split into two growth forms: trees and shrubs
divers_tspest <- summarise(
group_by(pestdat.ag_adult, site, plot, quadrat, quad.unique, growth.form, pesticide, exclosure, census),
shannon=diversity(survs), simpson=diversity(survs, index='invsimpson'), species=specnumber(survs))
divers_tspest$A.sum <- adult.ngbr.inf$A.sum[match(divers_tspest$quad.unique, adult.ngbr.inf$quad.unique)]
summary(divers_tspest)## site plot quadrat quad.unique growth.form pesticide
## 1:580 1_0:293 Min. :101.0 1_0_101: 8 shrub:753 W:400
## 2:572 1_1:287 1st Qu.:202.0 1_0_103: 8 tree :853 F:407
## 3:454 2_0:300 Median :304.0 1_0_105: 8 I:394
## 2_1:272 Mean :302.9 1_0_107: 8 C:405
## 3_0:229 3rd Qu.:408.0 1_0_108: 8
## 3_1:225 Max. :510.0 1_0_201: 8
## (Other):1558
## exclosure census shannon simpson species
## 0:822 16fa:426 Min. :0.0000 Min. :1.000 Min. :0.000
## 1:784 16sp:327 1st Qu.:0.0000 1st Qu.:1.324 1st Qu.:1.000
## 17fa:428 Median :0.6931 Median :2.000 Median :2.000
## 17sp:425 Mean :0.6933 Mean : Inf Mean :2.504
## 3rd Qu.:1.0735 3rd Qu.:2.830 3rd Qu.:3.000
## Max. :1.9560 Max. : Inf Max. :9.000
##
## A.sum
## Min. : 3990
## 1st Qu.: 4991
## Median : 5455
## Mean : 5630
## 3rd Qu.: 6045
## Max. :10411
## divers_tspest$A.sum_curt <- divers_tspest$A.sum^(1/3)
### Trees
# shannon diversity
mod.tree.shannon1 <- lmer(shannon ~ site + census + pesticide*exclosure + (1|quad.unique),
data=subset(divers_tspest, growth.form=='tree'))
mod.tree.shannon2 <- lmer(shannon ~ site + census + pesticide + exclosure + (1|quad.unique),
data=subset(divers_tspest, growth.form=='tree'))
mod.tree.shannon3 <- lmer(shannon ~ site + census + pesticide + exclosure + A.sum_curt +
(1|quad.unique), data=subset(divers_tspest, growth.form=='tree'))
anova(mod.tree.shannon1, mod.tree.shannon2, mod.tree.shannon3)## refitting model(s) with ML (instead of REML)## Data: subset(divers_tspest, growth.form == "tree")
## Models:
## mod.tree.shannon2: shannon ~ site + census + pesticide + exclosure + (1 | quad.unique)
## mod.tree.shannon3: shannon ~ site + census + pesticide + exclosure + A.sum_curt +
## mod.tree.shannon3: (1 | quad.unique)
## mod.tree.shannon1: shannon ~ site + census + pesticide * exclosure + (1 | quad.unique)
## Df AIC BIC logLik deviance Chisq Chi Df
## mod.tree.shannon2 12 516.04 573.03 -246.02 492.04
## mod.tree.shannon3 13 517.96 579.69 -245.98 491.96 0.0892 1
## mod.tree.shannon1 15 521.59 592.83 -245.80 491.59 0.3612 2
## Pr(>Chisq)
## mod.tree.shannon2
## mod.tree.shannon3 0.7652
## mod.tree.shannon1 0.8348plot(mod.tree.shannon2, type=c('p', 'smooth'))
qqPlot(resid(mod.tree.shannon2))
sjp.lmer(mod.tree.shannon2, type='fe')## Computing p-values via Kenward-Roger approximation. Use `p.kr = FALSE` if computation takes too long.
sjp.lmer(mod.tree.shannon2, type='re', sort.est="sort.all" )## Sorting each group of random intercept ('sort.all') is not possible when 'facet.grid = TRUE'.
## Plotting random effects...
summary(mod.tree.shannon2, correlation=FALSE)## Linear mixed model fit by REML ['lmerMod']
## Formula:
## shannon ~ site + census + pesticide + exclosure + (1 | quad.unique)
## Data: subset(divers_tspest, growth.form == "tree")
##
## REML criterion at convergence: 543.5
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.4107 -0.5661 0.0429 0.6062 2.9724
##
## Random effects:
## Groups Name Variance Std.Dev.
## quad.unique (Intercept) 0.06964 0.2639
## Residual 0.06892 0.2625
## Number of obs: 853, groups: quad.unique, 239
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 0.81838 0.04553 17.975
## site1 0.09859 0.02740 3.598
## site2 0.03608 0.02724 1.324
## census16sp -0.69200 0.02888 -23.958
## census17fa 0.23466 0.02410 9.735
## census17sp -0.10985 0.02412 -4.555
## pesticideF 0.12179 0.05440 2.239
## pesticideI 0.03059 0.05487 0.558
## pesticideC -0.06755 0.05471 -1.235
## exclosure1 -0.07763 0.03870 -2.006## Invsimpson
divers_tspest1 <- subset(divers_tspest, simpson != Inf)
mod.tree.simpson <- lmer(simpson ~ site + census + pesticide + exclosure + (1|quad.unique),
data=subset(divers_tspest1, growth.form=='tree'))
plot(mod.tree.simpson, type=c('p', 'smooth'))
qqPlot(resid(mod.tree.simpson))
sjp.lmer(mod.tree.simpson, type='fe')## Computing p-values via Kenward-Roger approximation. Use `p.kr = FALSE` if computation takes too long.
sjp.lmer(mod.tree.simpson, type='re', sort.est="sort.all" )## Sorting each group of random intercept ('sort.all') is not possible when 'facet.grid = TRUE'.
## Plotting random effects...
summary(mod.tree.simpson, correlation=FALSE)## Linear mixed model fit by REML ['lmerMod']
## Formula:
## simpson ~ site + census + pesticide + exclosure + (1 | quad.unique)
## Data: subset(divers_tspest1, growth.form == "tree")
##
## REML criterion at convergence: 1764.7
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.8480 -0.5775 -0.0833 0.5437 3.5601
##
## Random effects:
## Groups Name Variance Std.Dev.
## quad.unique (Intercept) 0.3053 0.5525
## Residual 0.3258 0.5707
## Number of obs: 815, groups: quad.unique, 238
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 2.16286 0.09681 22.340
## site1 0.23976 0.05823 4.118
## site2 0.02483 0.05789 0.429
## census16sp -0.99589 0.06912 -14.407
## census17fa 0.59961 0.05263 11.394
## census17sp -0.14007 0.05265 -2.660
## pesticideF 0.24869 0.11558 2.152
## pesticideI 0.12667 0.11689 1.084
## pesticideC -0.06748 0.11609 -0.581
## exclosure1 -0.21902 0.08230 -2.661#### shrub
# shannon diversity
mod.shrub.shannon1 <- lmer(shannon ~ site + census + pesticide*exclosure + (1|quad.unique),
data=subset(divers_tspest, growth.form=='shrub'))
mod.shrub.shannon2 <- lmer(shannon ~ site + census + pesticide + exclosure + (1|quad.unique),
data=subset(divers_tspest, growth.form=='shrub'))
mod.shrub.shannon3 <- lmer(shannon ~ site + census + pesticide + exclosure + A.sum_curt +
(1|quad.unique), data=subset(divers_tspest, growth.form=='shrub'))
anova(mod.shrub.shannon1, mod.shrub.shannon2, mod.shrub.shannon3)## refitting model(s) with ML (instead of REML)## Data: subset(divers_tspest, growth.form == "shrub")
## Models:
## mod.shrub.shannon2: shannon ~ site + census + pesticide + exclosure + (1 | quad.unique)
## mod.shrub.shannon3: shannon ~ site + census + pesticide + exclosure + A.sum_curt +
## mod.shrub.shannon3: (1 | quad.unique)
## mod.shrub.shannon1: shannon ~ site + census + pesticide * exclosure + (1 | quad.unique)
## Df AIC BIC logLik deviance Chisq Chi Df
## mod.shrub.shannon2 12 -191.34 -135.85 107.67 -215.34
## mod.shrub.shannon3 13 -191.40 -131.29 108.70 -217.40 2.0587 1
## mod.shrub.shannon1 15 -186.12 -116.76 108.06 -216.12 0.0000 2
## Pr(>Chisq)
## mod.shrub.shannon2
## mod.shrub.shannon3 0.1513
## mod.shrub.shannon1 1.0000plot(mod.shrub.shannon2, type=c('p', 'smooth'))
qqPlot(resid(mod.shrub.shannon2))
sjp.lmer(mod.shrub.shannon2, type='fe')## Computing p-values via Kenward-Roger approximation. Use `p.kr = FALSE` if computation takes too long.
sjp.lmer(mod.shrub.shannon2, type='re', sort.est="sort.all" )## Sorting each group of random intercept ('sort.all') is not possible when 'facet.grid = TRUE'.
## Plotting random effects...
summary(mod.shrub.shannon2, correlation=FALSE)## Linear mixed model fit by REML ['lmerMod']
## Formula:
## shannon ~ site + census + pesticide + exclosure + (1 | quad.unique)
## Data: subset(divers_tspest, growth.form == "shrub")
##
## REML criterion at convergence: -167.7
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -5.4242 -0.2175 -0.0030 0.1953 3.9673
##
## Random effects:
## Groups Name Variance Std.Dev.
## quad.unique (Intercept) 0.22503 0.4744
## Residual 0.01582 0.1258
## Number of obs: 753, groups: quad.unique, 195
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 0.6630352 0.0778860 8.513
## site1 0.1674967 0.0470239 3.562
## site2 0.0603374 0.0484612 1.245
## census16sp -0.0247005 0.0130640 -1.891
## census17fa 0.0286765 0.0130150 2.203
## census17sp 0.0002409 0.0129706 0.019
## pesticideF -0.0187711 0.0984434 -0.191
## pesticideI -0.0345274 0.0973471 -0.355
## pesticideC -0.0640912 0.0963271 -0.665
## exclosure1 -0.0697174 0.0688905 -1.012## Invsimpson
mod.shrub.simpson <- lmer(simpson ~ site + census + pesticide + exclosure + (1|quad.unique),
data=subset(divers_tspest1, growth.form=='shrub'))
plot(mod.shrub.simpson, type=c('p', 'smooth'))
qqPlot(resid(mod.shrub.simpson))
sjp.lmer(mod.shrub.simpson, type='fe')## Computing p-values via Kenward-Roger approximation. Use `p.kr = FALSE` if computation takes too long.
sjp.lmer(mod.shrub.simpson, type='re', sort.est="sort.all" )## Sorting each group of random intercept ('sort.all') is not possible when 'facet.grid = TRUE'.
## Plotting random effects...
summary(mod.shrub.simpson, correlation=FALSE)## Linear mixed model fit by REML ['lmerMod']
## Formula:
## simpson ~ site + census + pesticide + exclosure + (1 | quad.unique)
## Data: subset(divers_tspest1, growth.form == "shrub")
##
## REML criterion at convergence: 867.9
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.7200 -0.2548 -0.0147 0.1818 3.7375
##
## Random effects:
## Groups Name Variance Std.Dev.
## quad.unique (Intercept) 0.88945 0.9431
## Residual 0.06534 0.2556
## Number of obs: 747, groups: quad.unique, 192
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 2.114774 0.155246 13.622
## site1 0.320558 0.094192 3.403
## site2 0.053499 0.097229 0.550
## census16sp -0.051347 0.026575 -1.932
## census17fa 0.061410 0.026476 2.319
## census17sp -0.002876 0.026533 -0.108
## pesticideF -0.066060 0.198170 -0.333
## pesticideI -0.054336 0.194706 -0.279
## pesticideC -0.121789 0.191635 -0.636
## exclosure1 -0.092261 0.138194 -0.668When trees and shrubs are pooled in the diversity test, pesticide had no effect on diversity. Exclosure lower diversity. When trees and shrubs are estimated seperately, fungicide higher diversity of tree seedlings.